{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 对数据做数据探索分析（可参考EDA_BikeSharing.ipynb，不计分） \n",
    "2. 适当的特征工程（可参考FE_BikeSharing.ipynb，不计分） \n",
    "3. 对全体数据，随机选择其中80%做训练数据，剩下20%为测试数据，评价指标为RMSE。（10分） \n",
    "4. 用训练数据训练最小二乘线性回归模型（20分）、岭回归模型、Lasso模型，其中岭回归模型（30分）和Lasso模型（30分），注意岭回归模型和Lasso模型的正则超参数调优。 \n",
    "5. 比较用上述三种模型得到的各特征的系数，以及各模型在测试集上的性能。并简单说明原因。（10分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import r2_score, mean_squared_error  #评价回归预测模型的性能\n",
    "from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV\n",
    "from math import sqrt\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sn\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(584, 33)\n",
      "<bound method NDFrame.head of      season_1  season_2  season_3  season_4  mnth_1  mnth_2  mnth_3  mnth_4  \\\n",
      "375         1         0         0         0       1       0       0       0   \n",
      "32          1         0         0         0       0       1       0       0   \n",
      "591         0         0         1         0       0       0       0       0   \n",
      "376         1         0         0         0       1       0       0       0   \n",
      "676         0         0         0         1       0       0       0       0   \n",
      "328         0         0         0         1       0       0       0       0   \n",
      "381         1         0         0         0       1       0       0       0   \n",
      "273         0         0         0         1       0       0       0       0   \n",
      "116         0         1         0         0       0       0       0       1   \n",
      "317         0         0         0         1       0       0       0       0   \n",
      "10          1         0         0         0       1       0       0       0   \n",
      "3           1         0         0         0       1       0       0       0   \n",
      "530         0         1         0         0       0       0       0       0   \n",
      "178         0         0         1         0       0       0       0       0   \n",
      "13          1         0         0         0       1       0       0       0   \n",
      "257         0         0         1         0       0       0       0       0   \n",
      "416         1         0         0         0       0       1       0       0   \n",
      "40          1         0         0         0       0       1       0       0   \n",
      "226         0         0         1         0       0       0       0       0   \n",
      "65          1         0         0         0       0       0       1       0   \n",
      "72          1         0         0         0       0       0       1       0   \n",
      "284         0         0         0         1       0       0       0       0   \n",
      "650         0         0         0         1       0       0       0       0   \n",
      "113         0         1         0         0       0       0       0       1   \n",
      "707         0         0         0         1       0       0       0       0   \n",
      "512         0         1         0         0       0       0       0       0   \n",
      "335         0         0         0         1       0       0       0       0   \n",
      "182         0         0         1         0       0       0       0       0   \n",
      "623         0         0         1         0       0       0       0       0   \n",
      "399         1         0         0         0       0       1       0       0   \n",
      "..        ...       ...       ...       ...     ...     ...     ...     ...   \n",
      "325         0         0         0         1       0       0       0       0   \n",
      "323         0         0         0         1       0       0       0       0   \n",
      "255         0         0         1         0       0       0       0       0   \n",
      "570         0         0         1         0       0       0       0       0   \n",
      "545         0         0         1         0       0       0       0       0   \n",
      "550         0         0         1         0       0       0       0       0   \n",
      "359         1         0         0         0       0       0       0       0   \n",
      "51          1         0         0         0       0       1       0       0   \n",
      "102         0         1         0         0       0       0       0       1   \n",
      "195         0         0         1         0       0       0       0       0   \n",
      "679         0         0         0         1       0       0       0       0   \n",
      "628         0         0         1         0       0       0       0       0   \n",
      "124         0         1         0         0       0       0       0       0   \n",
      "404         1         0         0         0       0       1       0       0   \n",
      "287         0         0         0         1       0       0       0       0   \n",
      "47          1         0         0         0       0       1       0       0   \n",
      "587         0         0         1         0       0       0       0       0   \n",
      "607         0         0         1         0       0       0       0       0   \n",
      "433         1         0         0         0       0       0       1       0   \n",
      "674         0         0         0         1       0       0       0       0   \n",
      "263         0         0         1         0       0       0       0       0   \n",
      "360         1         0         0         0       0       0       0       0   \n",
      "75          1         0         0         0       0       0       1       0   \n",
      "466         0         1         0         0       0       0       0       1   \n",
      "299         0         0         0         1       0       0       0       0   \n",
      "534         0         1         0         0       0       0       0       0   \n",
      "584         0         0         1         0       0       0       0       0   \n",
      "493         0         1         0         0       0       0       0       0   \n",
      "527         0         1         0         0       0       0       0       0   \n",
      "168         0         1         0         0       0       0       0       0   \n",
      "\n",
      "     mnth_5  mnth_6  ...  weekday_4  weekday_5  weekday_6      temp     atemp  \\\n",
      "375       0       0  ...          0          0          0  0.267946  0.267451   \n",
      "32        0       0  ...          0          0          0  0.250293  0.230167   \n",
      "591       0       0  ...          0          0          0  0.831783  0.784574   \n",
      "376       0       0  ...          1          0          0  0.402934  0.397556   \n",
      "676       0       0  ...          0          0          0  0.294943  0.256688   \n",
      "328       0       0  ...          0          1          0  0.393589  0.395893   \n",
      "381       0       0  ...          0          0          0  0.391151  0.374375   \n",
      "273       0       0  ...          0          0          1  0.437201  0.438149   \n",
      "116       0       0  ...          0          0          0  0.698871  0.651162   \n",
      "317       0       0  ...          0          0          0  0.586727  0.584835   \n",
      "10        0       0  ...          0          0          0  0.137016  0.147533   \n",
      "3         0       0  ...          0          0          0  0.175530  0.174649   \n",
      "530       0       1  ...          1          0          0  0.734175  0.715797   \n",
      "178       0       1  ...          0          0          0  0.853589  0.805286   \n",
      "13        0       0  ...          0          1          0  0.126773  0.143528   \n",
      "257       0       0  ...          1          0          0  0.645914  0.622978   \n",
      "416       0       0  ...          0          0          0  0.284966  0.283586   \n",
      "40        0       0  ...          1          0          0  0.106185  0.092512   \n",
      "226       0       0  ...          0          0          0  0.755981  0.705013   \n",
      "65        0       0  ...          0          0          0  0.252460  0.209223   \n",
      "72        0       0  ...          0          0          0  0.331557  0.332005   \n",
      "284       0       0  ...          0          0          0  0.603340  0.575784   \n",
      "650       0       0  ...          0          1          0  0.471467  0.464675   \n",
      "113       0       0  ...          0          0          0  0.651106  0.620474   \n",
      "707       0       0  ...          0          0          1  0.401896  0.407492   \n",
      "512       1       0  ...          0          0          0  0.786094  0.738167   \n",
      "335       0       0  ...          0          1          0  0.317788  0.331261   \n",
      "182       0       0  ...          0          0          1  0.846320  0.772142   \n",
      "623       0       0  ...          0          0          1  0.684333  0.665240   \n",
      "399       0       0  ...          0          0          1  0.255486  0.254199   \n",
      "..      ...     ...  ...        ...        ...        ...       ...       ...   \n",
      "325       0       0  ...          0          0          0  0.445508  0.449743   \n",
      "323       0       0  ...          0          0          0  0.503656  0.496161   \n",
      "255       0       0  ...          0          0          0  0.737290  0.688457   \n",
      "570       0       0  ...          0          0          0  0.861895  0.823521   \n",
      "545       0       1  ...          0          1          0  0.965734  0.928746   \n",
      "550       0       0  ...          0          0          0  0.909661  0.857502   \n",
      "359       0       0  ...          0          0          0  0.327223  0.310393   \n",
      "51        0       0  ...          0          0          0  0.304288  0.269097   \n",
      "102       0       0  ...          0          0          0  0.440316  0.443951   \n",
      "195       0       0  ...          0          1          0  0.752866  0.715782   \n",
      "679       0       0  ...          0          0          1  0.411242  0.413306   \n",
      "628       0       0  ...          1          0          0  0.607495  0.594784   \n",
      "124       1       0  ...          1          0          0  0.498465  0.476286   \n",
      "404       0       0  ...          1          0          0  0.256524  0.240105   \n",
      "287       0       0  ...          0          0          1  0.557653  0.550854   \n",
      "47        0       0  ...          1          0          0  0.469390  0.458882   \n",
      "587       0       0  ...          0          1          0  0.818284  0.772975   \n",
      "607       0       0  ...          1          0          0  0.806862  0.753071   \n",
      "433       0       0  ...          0          1          0  0.438239  0.417436   \n",
      "674       0       0  ...          0          0          0  0.324018  0.300601   \n",
      "263       0       0  ...          0          0          0  0.667720  0.618859   \n",
      "360       0       0  ...          0          0          0  0.331287  0.326273   \n",
      "75        0       0  ...          1          0          0  0.443431  0.434828   \n",
      "466       0       0  ...          0          0          0  0.360813  0.338928   \n",
      "299       0       0  ...          1          0          0  0.511964  0.496145   \n",
      "534       0       1  ...          0          0          0  0.634491  0.611389   \n",
      "584       0       0  ...          0          0          0  0.843205  0.811932   \n",
      "493       1       0  ...          0          0          0  0.651106  0.627966   \n",
      "527       0       1  ...          0          0          0  0.824514  0.762183   \n",
      "168       0       1  ...          0          0          1  0.794402  0.741487   \n",
      "\n",
      "          hum  windspeed  holiday  workingday  yr  \n",
      "375  0.871465   0.224357        0           1   1  \n",
      "32   0.797344   0.498723        0           1   0  \n",
      "591  0.706084   0.302566        0           1   1  \n",
      "376  0.825622   0.326911        0           1   1  \n",
      "676  0.562982   0.580773        0           1   1  \n",
      "328  0.661954   0.157717        0           1   0  \n",
      "381  0.736336   0.673588        0           1   1  \n",
      "273  0.775064   0.556422        0           0   0  \n",
      "116  0.859041   0.597455        0           1   0  \n",
      "317  0.603684   0.585902        0           1   0  \n",
      "10   0.705773   0.205620        0           1   0  \n",
      "3    0.607131   0.284297        0           1   0  \n",
      "530  0.585689   0.476922        0           1   1  \n",
      "178  0.652100   0.251285        0           1   0  \n",
      "13   0.553034   0.214724        0           1   0  \n",
      "257  0.729221   0.512820        0           1   0  \n",
      "416  0.611155   0.377935        0           1   1  \n",
      "40   0.449759   0.411369        0           1   0  \n",
      "226  0.732219   0.384608        0           1   0  \n",
      "65   0.566894   0.657553        0           1   0  \n",
      "72   0.511010   0.236118        0           1   0  \n",
      "284  0.931877   0.465413        0           1   0  \n",
      "650  0.554413   0.438493        0           1   1  \n",
      "113  0.833761   0.350017        0           0   0  \n",
      "707  0.937018   0.162836        0           0   1  \n",
      "512  0.716795   0.397425        0           0   1  \n",
      "335  0.643530   0.161548        0           1   0  \n",
      "182  0.457155   0.191045        0           0   0  \n",
      "623  0.515853   0.464116        0           0   1  \n",
      "399  0.801628   0.205133        0           0   1  \n",
      "..        ...        ...      ...         ...  ..  \n",
      "325  0.989717   0.198734        0           1   0  \n",
      "323  0.703941   0.337184        0           0   0  \n",
      "255  0.732648   0.246175        0           1   0  \n",
      "570  0.673522   0.389762        0           1   1  \n",
      "545  0.502571   0.294854        0           1   1  \n",
      "550  0.553985   0.232045        1           0   1  \n",
      "359  0.521293   0.447508        1           0   0  \n",
      "51   0.622108   0.588479        1           0   0  \n",
      "102  0.842331   0.470498        0           1   0  \n",
      "195  0.607969   0.330758        0           1   0  \n",
      "679  0.663668   0.073090        0           0   1  \n",
      "628  0.635818   0.197446        0           1   1  \n",
      "124  0.456727   0.562804        0           1   0  \n",
      "404  0.577977   0.353856        0           1   1  \n",
      "287  0.497001   0.485890        0           0   0  \n",
      "47   0.519280   0.428210        0           1   0  \n",
      "587  0.736075   0.446164        0           1   1  \n",
      "607  0.607113   0.112818        0           1   1  \n",
      "433  0.418594   0.808970        0           1   1  \n",
      "674  0.508141   0.441035        0           1   1  \n",
      "263  0.925450   0.152581        0           1   0  \n",
      "360  0.784062   0.342338        0           1   0  \n",
      "75   0.619966   0.385896        0           1   0  \n",
      "466  0.482843   0.562561        0           1   1  \n",
      "299  0.835904   0.361537        0           1   0  \n",
      "534  0.799915   0.314086        0           1   1  \n",
      "584  0.723650   0.194850        0           1   1  \n",
      "493  0.705227   0.564118        0           1   1  \n",
      "527  0.604542   0.382050        0           1   1  \n",
      "168  0.689375   0.200004        0           0   0  \n",
      "\n",
      "[584 rows x 33 columns]>\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 读取去量纲后的数据\n",
    "df = pd.read_csv(\"bike_day_FE.csv\")\n",
    "y = df[\"cnt\"]\n",
    "X = df.drop([\"cnt\", \"instant\"], axis=1)\n",
    "feat_names = X.columns\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y,random_state=2,test_size=0.2)\n",
    "print(X_train.shape)\n",
    "print(X_train.head)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def RMSE(targets, preds):\n",
    "    return sqrt(mean_squared_error( targets , preds))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         columns          coef\n",
      "16  weathersit_1  1.749291e+16\n",
      "17  weathersit_2  1.749291e+16\n",
      "18  weathersit_3  1.749291e+16\n",
      "25     weekday_6  2.402360e+15\n",
      "19     weekday_0  2.402360e+15\n",
      "31    workingday  1.998868e+15\n",
      "30       holiday  1.998868e+15\n",
      "24     weekday_5  4.034916e+14\n",
      "23     weekday_4  4.034916e+14\n",
      "22     weekday_3  4.034916e+14\n",
      "21     weekday_2  4.034916e+14\n",
      "20     weekday_1  4.034916e+14\n",
      "26          temp  2.785642e+03\n",
      "32            yr  1.977924e+03\n",
      "27         atemp  1.093423e+03\n",
      "3       season_4  7.199798e+02\n",
      "1       season_2  6.450353e+01\n",
      "2       season_3  4.123943e+01\n",
      "0       season_1 -8.709105e+02\n",
      "28           hum -1.518035e+03\n",
      "29     windspeed -1.534058e+03\n",
      "12        mnth_9 -1.335143e+16\n",
      "8         mnth_5 -1.335143e+16\n",
      "6         mnth_3 -1.335143e+16\n",
      "7         mnth_4 -1.335143e+16\n",
      "13       mnth_10 -1.335143e+16\n",
      "9         mnth_6 -1.335143e+16\n",
      "11        mnth_8 -1.335143e+16\n",
      "5         mnth_2 -1.335143e+16\n",
      "4         mnth_1 -1.335143e+16\n",
      "15       mnth_12 -1.335143e+16\n",
      "14       mnth_11 -1.335143e+16\n",
      "10        mnth_7 -1.335143e+16\n",
      "The RMSE of LinearRegression on test is 798.511536\n",
      "The RMSE of LinearRegression on train is 745.108925\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAFsCAYAAADyj6FyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAeM0lEQVR4nO3df7hWZZ3v8fdX2EKJoCgYQgWWVqgJtjWZ0igqLUyc61TWMcWRtKPVsTpW/phGxmquGp2aPDk1nH6oaYkyUzZ1NelhZExPaWBkKhqkUiADqJlimYDf88dawGazN/vZv9g3+3m/ruu5nrXu9ete94b9ee617r2eyEwkSdLA22OgKyBJkiqGsiRJhTCUJUkqhKEsSVIhDGVJkgphKEuSVAhDWU0vIu6LiOkDXY+BFBF/GRG/i4gNETF1Fx53Q0Qc1M/HyIh4eX8eQ+orhrIGtYh4JCLe3K7sjIi4fct8Zh6amYu62M/E+pf70H6q6kC7HPhQZo7IzF+0X1if+zN1iK6OiC9ExJDeHrQ+3kO93Y80WBjKUgEKCPuXAvd1sc4RmTkCeANwCnBmv9dKajKGsppe2950RBwdEYsj4qmIWBsRX6hXu61+f7LuLU6LiD0i4q8jYmVErIuIayJiVJv9nl4vezwiPtXuOHMjYkFEXBsRTwFn1Mf+aUQ8GRFrIuLLEbFnm/1lRJwbEcsj4umI+HREvKze5qmIuKHt+u3OscO6RsSwiNgADAF+GRG/6aq9MnMFcAcwpc3+R0XE1+t6r46Iz2zpSUfEyyPiPyPiDxHxWETMb3dOL6+nXxAR/1DX8Q8RcXtdNj0iVnXxM+u03dpt9/aIuL9uv9URcX5X5yvtSoaytL0vAV/KzJHAy4Ab6vLj6vd96kuuPwXOqF9vBA4CRgBfBoiIycA/AacC44BRwPh2x5oFLAD2Aa4DNgMfBfYHpgEzgHPbbXMC8BrgGOATwLz6GC8GDgPe28l5dVjXzPxz3fuFqif8ss6bphIRrwSOBVa0Kb4a2AS8HJgKvBV4f73s08DNwL7ABOB/d7Lry+tz+wtgdH1+z3dVHxprty2+DnwgM/emaq//aGD/0i5jKKsZfK/uRT0ZEU9ShWVnNgIvj4j9M3NDZv5sJ+ueCnwhMx/KzA3AhcB76kvR7wT+LTNvz8zngL8B2j9o/qeZ+b3MfD4z/5SZSzLzZ5m5KTMfAf6Z6lJxW5/PzKcy8z7gXuDm+vh/AH5EFYjdrWuj7o6IZ4BlwCLqdoyIA4C3AR/JzGcycx3wReA99XYbqS6PH5iZz2bm7e13HBF7UF0OPy8zV2fm5sz8f5n5564q1WC7bbERmBwRIzPz95l5d+OnL/U/Q1nN4OTM3GfLi857UQBzgEOAByLi5xFx4k7WPRBY2WZ+JTAUOKBe9rstCzLzj8Dj7bb/XduZiDgkIn4QEf9VX9L+O6reX1tr20z/qYP5EXRsZ3Vt1JH1/k8BXgvsVZe/FGgB1rT54PPPwNh6+SeAAO6KaqR7R/ei9weGA11ePm+vwXbb4r8BbwdW1pfUp3X3eFJ/MpSlNjJzeWa+lypQPg8siIi92LGXC/AoVSBt8RKqS7hrgTVUl2qB6n4psF/7w7Wb/wrwAHBwffn8Iqow6ws7q2vDsnID8FOq3j9UHy7+DOzf5sPPyMw8tN7mvzLzrMw8EPgA8E+x458oPQY8S3XLoL1ngBdumanvVY9ps7zhdsvMn2fmLKqf7/fYdntCKoKhLLUREe+LiDGZ+TzwZF28GVhPdX+z7d/Ufgf4aERMiogRVD20+Zm5iepe8Tsi4i/qQUd/S9cBuzfwFLChvm97Tp+d2M7r2hOfA86OiBdl5hqqe8b/EBEj60FlL4uINwBExLsiYssHlN9TfRjZ3HZndXt/A/hCRBwYEUOiGkw3DPg1MDwiZkZEC/DXwLA2mzfUbhGxZ0ScGhGjMnNjvc3mjtaVBoqhLG3vBOC+ekTyl4D31PdB/wh8FrijvkR7DFWIfItqZPbDVD29DwPU93w/DFxP1Wt+GlhH1aPszPnAf6/X/T/A/J2s212d1rUnMvNXwH8CH6+LTgf2BO6nCt4FVAPcAI4C7qzb9PtU940f7mC35wO/An4OPEF1pWKP+n75ucDXgNVUPedV7bZrtN1OAx6pL3P/D+B9jZ+11P8is6OrcpL6Ut07fZLqEmtHgSRJ9pSl/hIR74iIF9b3pC+n6gU+MrC1klQyQ1nqP7OoBlg9ChxMdSncS1OSOuXla0mSCmFPWZKkQuzSh+Dvv//+OXHixF15SEmSirJkyZLHMnNMR8t2aShPnDiRxYsX78pDSpJUlIhY2dkyL19LklQIQ1mSpEIYypIkFWKX3lOWJHVs48aNrFq1imeffXagq6I+Mnz4cCZMmEBLS0vD2xjKklSAVatWsffeezNx4kQi+urLwTRQMpPHH3+cVatWMWnSpIa38/K1JBXg2WefZb/99jOQB4mIYL/99uv2lQ9DWZIKYSAPLj35eRrKkiQVwnvKklSguXN3/f6GDBnC4YcfzqZNm5g0aRLf+ta32Geffbp9rPe///187GMfY/LkyduVX3XVVSxevJgvf/nL3d4nwIgRI9iwYUND606fPp3LL7+c1tbWrWWLFy/mmmuu4YorrujR8XcFe8qSJABe8IIXsHTpUu69915Gjx7NlVde2aP9fO1rX9shkEvQ2tra74G8efPmXm1vKEuSdjBt2jRWr169df6yyy7jqKOO4tWvfjWXXHIJAM888wwzZ87kiCOO4LDDDmP+/PlA1Uvd8kjlb37zmxxyyCG84Q1v4I477ti6vzPOOIMFCxZsnR8xYgQAGzZsYMaMGRx55JEcfvjh3HTTTTvUbc2aNRx33HFMmTKFww47jJ/85CcNndOiRYs48cQTAZg7dy5nnnkm06dP56CDDtourK+99lqOPvpopkyZwgc+8IGtQXvOOefQ2trKoYceurUNoHqE9KWXXsrrX/96brzxxobq0hkvX0uStrN582YWLlzInDlzALj55ptZvnw5d911F5nJSSedxG233cb69es58MAD+eEPfwjAH/7wh+32s2bNGi655BKWLFnCqFGjeOMb38jUqVN3euzhw4fz3e9+l5EjR/LYY49xzDHHcNJJJ203aOrb3/42xx9/PBdffDGbN2/mj3/8Y4/O84EHHuDWW2/l6aef5hWveAXnnHMOK1asYP78+dxxxx20tLRw7rnnct1113H66afz2c9+ltGjR7N582ZmzJjBPffcw6tf/eqt9b799tt7VI+2DGVJEgB/+tOfmDJlCo888givec1reMtb3gJUoXzzzTdvDdQNGzawfPlyjj32WM4//3w++clPcuKJJ3Lsscdut78777yT6dOnM2ZM9YVIp5xyCr/+9a93WofM5KKLLuK2225jjz32YPXq1axdu5YXvehFW9c56qijOPPMM9m4cSMnn3wyU6ZM6dH5zpw5k2HDhjFs2DDGjh3L2rVrWbhwIUuWLOGoo47a2iZjx44F4IYbbmDevHls2rSJNWvWcP/9928N5VNOOaVHdWjPy9eSJGDbPeWVK1fy3HPPbb2nnJlceOGFLF26lKVLl7JixQrmzJnDIYccwpIlSzj88MO58MILufTSS3fYZ2d/FjR06FCef/75rft/7rnnALjuuutYv349S5YsYenSpRxwwAE7/K3vcccdx2233cb48eM57bTTuOaaa3p0vsOGDds6PWTIEDZt2kRmMnv27K3n+uCDDzJ37lwefvhhLr/8chYuXMg999zDzJkzt6vXXnvt1aM6tGdPWdJuq9ERyn09knmwGzVqFFdccQWzZs3inHPO4fjjj+dTn/oUp556KiNGjGD16tW0tLSwadMmRo8ezfve9z5GjBjBVVddtd1+Xvva13Leeefx+OOPM3LkSG688UaOOOIIoLoPu2TJEt797ndz0003sXHjRqC6BD527FhaWlq49dZbWblyx285XLlyJePHj+ess87imWee4e677+b000/vk3OfMWMGs2bN4qMf/Shjx47liSee4Omnn+app55ir732YtSoUaxdu5Yf/ehHTJ8+vU+O2ZahLEkFGugPElOnTuWII47g+uuv57TTTmPZsmVMmzYNqAZlXXvttaxYsYKPf/zj7LHHHrS0tPCVr3xlu32MGzeOuXPnMm3aNMaNG8eRRx65ddDUWWedxaxZszj66KOZMWPG1p7mqaeeyjve8Q5aW1uZMmUKr3zlK3eo26JFi7jssstoaWlhxIgRnfaUZ86cufW509OmTeODH/xgl+c9efJkPvOZz/DWt76V559/npaWFq688kqOOeYYpk6dyqGHHspBBx3E6173usYbsxsiM/tlxx1pbW3NLSPyJKm3BlNPedmyZbzqVa8a6Gqoj3X0c42IJZnZ2tH63lOWJKkQXr6WVJzdoWcr9Qd7ypJUiF15O1H9ryc/z4ZCOSL2iYgFEfFARCyLiGkRMToibomI5fX7vt0+uiQJqB4+8fjjjxvMg8SW71MePnx4t7Zr9PL1l4B/z8x3RsSewAuBi4CFmfm5iLgAuAD4ZLeOLkkCYMKECaxatYr169cPdFXUR4YPH86ECRO6tU2XoRwRI4HjgDMAMvM54LmImAVMr1e7GliEoSxJPdLS0sKkSZMGuhoaYI1cvj4IWA98MyJ+ERFfi4i9gAMycw1A/T62o40j4uyIWBwRi/0EKElS5xoJ5aHAkcBXMnMq8AzVpeqGZOa8zGzNzNYtzz+VJEk7aiSUVwGrMvPOen4BVUivjYhxAPX7uv6poiRJzaHLUM7M/wJ+FxGvqItmAPcD3wdm12WzgR2/9FKSJDWs0dHXHwauq0dePwT8FVWg3xARc4DfAu/qnypKktQcGgrlzFwKdPSczhl9Wx1JkpqXT/SSJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEI0+ZlOSemXu3IGugVQ+e8qSJBXCUJYkqRCGsiRJhTCUJUkqhKEsSVIhDGVJkgphKEuSVAhDWZKkQhjKkiQVwid6SRr0Gn2amE8d00CzpyxJUiEMZUmSCmEoS5JUCENZkqRCONBLkmrdGejloDD1B3vKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUiIa+TzkiHgGeBjYDmzKzNSJGA/OBicAjwLsz8/f9U01Jkga/7vSU35iZUzKztZ6/AFiYmQcDC+t5SZLUQ725fD0LuLqevho4uffVkSSpeTUaygncHBFLIuLsuuyAzFwDUL+P7WjDiDg7IhZHxOL169f3vsaSJA1SDd1TBl6XmY9GxFjgloh4oNEDZOY8YB5Aa2tr9qCOkiQ1hYZ6ypn5aP2+DvgucDSwNiLGAdTv6/qrkpIkNYMuQzki9oqIvbdMA28F7gW+D8yuV5sN3NRflZQkqRk0cvn6AOC7EbFl/W9n5r9HxM+BGyJiDvBb4F39V01Jkga/LkM5Mx8Cjuig/HFgRn9USpKkZuQTvSRJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUCENZkqRCGMqSJBXCUJYkqRCGsiRJhTCUJUkqhKEsSVIhDGVJkgphKEuSVAhDWZKkQhjKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUCENZkqRCGMqSJBXCUJYkqRCGsiRJhWg4lCNiSET8IiJ+UM9Piog7I2J5RMyPiD37r5qSJA1+3ekpnwcsazP/eeCLmXkw8HtgTl9WTJKkZtNQKEfEBGAm8LV6PoA3AQvqVa4GTu6PCkqS1CyGNrjePwKfAPau5/cDnszMTfX8KmB8RxtGxNnA2QAveclLel5TSUWaO3egayANHl32lCPiRGBdZi5pW9zBqtnR9pk5LzNbM7N1zJgxPaymJEmDXyM95dcBJ0XE24HhwEiqnvM+ETG07i1PAB7tv2pKkjT4ddlTzswLM3NCZk4E3gP8R2aeCtwKvLNebTZwU7/VUpKkJtCbv1P+JPCxiFhBdY/5631TJUmSmlOjA70AyMxFwKJ6+iHg6L6vkiRJzcknekmSVAhDWZKkQhjKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUCENZkqRCDB3oCkgq09y5A10DqfnYU5YkqRCGsiRJhTCUJUkqhKEsSVIhHOglST3Q1wPhHFgnsKcsSVIxDGVJkgphKEuSVAhDWZKkQhjKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmF6DKUI2J4RNwVEb+MiPsi4m/r8kkRcWdELI+I+RGxZ/9XV5KkwauRnvKfgTdl5hHAFOCEiDgG+Dzwxcw8GPg9MKf/qilJ0uDXZShnZUM921K/EngTsKAuvxo4uV9qKElSk2jonnJEDImIpcA64BbgN8CTmbmpXmUVML5/qihJUnNoKJQzc3NmTgEmAEcDr+potY62jYizI2JxRCxev359z2sqSdIg163R15n5JLAIOAbYJyKG1osmAI92ss28zGzNzNYxY8b0pq6SJA1qjYy+HhMR+9TTLwDeDCwDbgXeWa82G7ipvyopSVIzGNr1KowDro6IIVQhfkNm/iAi7geuj4jPAL8Avt6P9ZQkadDrMpQz8x5gagflD1HdX5YkSX3AJ3pJklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUCENZkqRCGMqSJBXCUJYkqRCGsiRJhTCUJUkqhKEsSVIhDGVJkgphKEuSVAhDWZKkQhjKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFGDrQFZAkwdy5/bOudi/2lCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmF8OEhUpPxwRNSuewpS5JUiC5DOSJeHBG3RsSyiLgvIs6ry0dHxC0Rsbx+37f/qytJ0uDVSE95E/C/MvNVwDHAByNiMnABsDAzDwYW1vOSJKmHugzlzFyTmXfX008Dy4DxwCzg6nq1q4GT+6uSkiQ1g27dU46IicBU4E7ggMxcA1VwA2M72ebsiFgcEYvXr1/fu9pKkjSINRzKETEC+BfgI5n5VKPbZea8zGzNzNYxY8b0pI6SJDWFhkI5IlqoAvm6zPzXunhtRIyrl48D1vVPFSVJag6NjL4O4OvAssz8QptF3wdm19OzgZv6vnqSJDWPRh4e8jrgNOBXEbG0LrsI+BxwQ0TMAX4LvKt/qihJUnPoMpQz83YgOlk8o2+rI0lS8/KJXpIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUCENZkqRCGMqSJBWikcdsSirc3LkDXQNJfcGesiRJhTCUJUkqhKEsSVIhDGVJkgrhQC+pYA7gkpqLPWVJkgphKEuSVAhDWZKkQhjKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSqEoSxJUiEMZUmSCmEoS5JUiKEDXQFpMJk7t2/Xk9Rc7ClLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSpEl6EcEd+IiHURcW+bstERcUtELK/f9+3fakqSNPg10lO+CjihXdkFwMLMPBhYWM9LkqRe6DKUM/M24Il2xbOAq+vpq4GT+7hekiQ1nZ7eUz4gM9cA1O9jO1sxIs6OiMURsXj9+vU9PJwkSYNfvw/0ysx5mdmama1jxozp78NJkrTb6mkor42IcQD1+7q+q5IkSc2pp6H8fWB2PT0buKlvqiNJUvNq5E+ivgP8FHhFRKyKiDnA54C3RMRy4C31vCRJ6oUuv7oxM9/byaIZfVwXSVID/IrQwcsnekmSVAhDWZKkQhjKkiQVwlCWJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEIYypIkFcJQliSpEIayJEmFMJQlSSpEl1/dKA1WA/n1d36lnnaFvv535r/b/mdPWZKkQhjKkiQVwlCWJKkQ3lOWJPW5gRyzsTuzpyxJUiEMZUmSCmEoS5JUCENZkqRCONBL6oIDUSTtKvaUJUkqhKEsSVIhDGVJkgphKEuSVAgHeqlf+A1MktR99pQlSSqEoSxJUiEMZUmSCmEoS5JUCAd6abfhAC5pYA3k/8Fm+SpIe8qSJBXCUJYkqRCGsiRJhTCUJUkqxG490Ks7N/T7epCAT6LqG814zpLKUOLgMXvKkiQVolehHBEnRMSDEbEiIi7oq0pJktSMehzKETEEuBJ4GzAZeG9ETO6rikmS1Gx601M+GliRmQ9l5nPA9cCsvqmWJEnNJzKzZxtGvBM4ITPfX8+fBrw2Mz/Ubr2zgbPr2VcAD+5kt/sDj/WoQoOL7bCNbVGxHSq2Q8V22GZ3bIuXZuaYjhb0ZvR1dFC2Q8Jn5jxgXkM7jFicma29qNOgYDtsY1tUbIeK7VCxHbYZbG3Rm8vXq4AXt5mfADzau+pIktS8ehPKPwcOjohJEbEn8B7g+31TLUmSmk+PL19n5qaI+BDwY2AI8I3MvK+X9WnoMncTsB22sS0qtkPFdqjYDtsMqrbo8UAvSZLUt3yilyRJhTCUJUkqxICEckScHxEZEfvX8xERV9SP67wnIo5ss+7siFhev2a3KX9NRPyq3uaKiOjoT7SKFBGfrs9zaUTcHBEH1uXN1g6XRcQD9bl+NyL2abPswvqcHoyI49uUd/ho13rA4Z11+8yvBx/uFiLiXRFxX0Q8HxGt7ZY1TTt0ZbA/1jcivhER6yLi3jZloyPilvrneUtE7FuXd/t3xe4iIl4cEbdGxLL6/8V5dXlztEVm7tIX1Z9R/RhYCexfl70d+BHV3z4fA9xZl48GHqrf962n962X3QVMq7f5EfC2XX0uvWiDkW2m/yfw1SZth7cCQ+vpzwOfr6cnA78EhgGTgN9QDSYcUk8fBOxZrzO53uYG4D319FeBcwb6/LrRDq+ierDOIqC1TXlTtUMXbdTpOQ+WF3AccCRwb5uyvwcuqKcvaPN/pNu/K3aXFzAOOLKe3hv4df1/oSnaYiB6yl8EPsH2DxqZBVyTlZ8B+0TEOOB44JbMfCIzfw/cApxQLxuZmT/NqvWvAU7etafRc5n5VJvZvdjWFs3WDjdn5qZ69mdUf+sOVTtcn5l/zsyHgRVUj3Xt8NGu9dWBNwEL6u2vZvdqh2WZ2dGT7pqqHbow6B/rm5m3AU+0K55F9XOE7X+e3fpd0f+17zuZuSYz766nnwaWAeNpkrbYpaEcEScBqzPzl+0WjQd+12Z+VV22s/JVHZTvNiLisxHxO+BU4G/q4qZrhzbOpPq0C91vh/2AJ9sE/O7cDm3ZDtt0ds6D3QGZuQaqsALG1uXd/bexW4qIicBU4E6apC1685jNDkXE/wVe1MGii4GLqC5Z7rBZB2XZg/Ji7KwdMvOmzLwYuDgiLgQ+BFxCE7ZDvc7FwCbgui2bdbB+0vGHyEHTDh1t1kHZbt0OvTCYz60ndtvfCY2KiBHAvwAfycyndjJcZlC1RZ+Hcma+uaPyiDic6r7YL+vGnQDcHRFH0/kjO1cB09uVL6rLJ3SwfjE6a4cOfBv4IVUoN1071IMvTgRm1JfgYeePcO2o/DGqS1ZD617ibtcOnRh07dALzfpY37URMS4z19SXZNfV5d39XbFbiYgWqkC+LjP/tS5ujrYYqJvZwCNsG+g1k+1v1N+V227UP0x1k37fenp0vezn9bpbBji9faDOpQfnfnCb6Q8DC5q0HU4A7gfGtCs/lO0HOD1ENdBnaD09iW2DfQ6tt7mR7Qc4nTvQ59eD9ljE9gO9mrIdOmmbTs95ML2AiWw/0Osyth/c9Pf1dLd/V+wur/qcrgH+sV15U7TFQDb8I2wL5QCupBpd+at2v5jOpBrgsgL4qzblrcC99TZfpn462e7wovoEeC9wD/BvwPgmbYcVVPd8ltavr7ZZdnF9Tg/SZkQ51UjLX9fLLm5TfhDVSPQVdTANG+jz60Y7/CXVp/o/A2uBHzdjOzTQTh2e82B5Ad8B1gAb638Pc6jGCSwEltfvWz6Md/t3xe7yAl5PdZn5nja/G97eLG3hYzYlSSqET/SSJKkQhrIkSYUwlCVJKoShLElSIQxlSZIKYShLklQIQ1mSpEL8f5hi8lldtFT+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAADQCAYAAADcQn7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dfZRcZZngf093Kkl3BLoD0U06ySQ4THJESAJRonh2TFAiotCHDxOWmYkfK+uoK7AaTdQzgLJjnMiA7ioOK3pQUMJHDAEcI0vC7gyzARI7IQaSIRpI0kEJJg2atKS6+9k/7nub29X3q6rr1kfX8zunT9d9671VT93u+9T7Pp+iqhiGYWRBU7UFMAxj9GIKxjCMzDAFYxhGZpiCMQwjM0zBGIaRGWOqLUAWnHLKKTpjxoxqi2EYDcPWrVtfVtVJheOjUsHMmDGDLVu2VFsMw2gYROSFsHHbIhmGkRmZKhgRuVZEdorIr0TkJyIyXkRmisgTIvKciKwRkbFu7jh3vMc9PyPwOivd+G4RWZylzIZhlI/MFIyIdACfAear6luBZmAp8HXgZlU9DTgCfMyd8jHgiKr+OXCzm4eIvMWddzrwPuA7ItKcldyGYZSPrLdIY4AWERkDtAIvAouA+9zzdwCd7vHF7hj3/HkiIm78blV9TVX3AnuAt2cstzEKWdfVzbmrNjJzxcOcu2oj67q6qy3SqCczBaOq3cA3gH14iuUVYCvQo6p9btoBoMM97gD2u3P73PyTg+Mh5wwiIleJyBYR2XLo0KHyfyCjrlnX1c3KtTvo7ulFge6eXlau3WFKJmOy3CK1460+ZgJTgAnABSFT/WxLiXguanzogOptqjpfVedPmjTMW2Y0OKs37KY33z9krDffz+oNu6skUWOQpZv6PcBeVT0EICJrgXcCbSIyxq1SpgIH3fwDwDTggNtSnQQcDoz7BM8xjFQc7OlNPb6uq5vVG3ZzsKeXKW0tLF88i855wxbNRgqytMHsAxaISKuzpZwHPANsAi5zc5YBD7jH690x7vmN6tWSWA8sdV6mmcBpwJMZym2MQqa0taQar6etVD3YlDJbwajqEyJyH/BLoA/oAm4DHgbuFpEb3djt7pTbgR+JyB68lctS9zo7ReQePOXUB3xKVYeudQ0jgeWLZ7Fy7Y4h26SWXDPLF88aMi9uK9U5r6NmVje+IvRl9RUhMEyeasoso7Hg1Pz589UieUcvpd4wac6bueLh4QY+PEPgzUvmhiqpr11yRiY3bJy8567aSHfI9q6jrYXHVywa8hqVkFlEtqrq/MLxUZkqYIxeivnmLqRzXkfinCltLaE37kktOT57z3b6C76Qg6ubNLKHKYywcWDY57x2zTa2vHCYGzvPCJURhtuUklZkWWMKxqgrsrxh1nV1c+x437DxXJNw9HjfMOXiE3WzF752mGLc8sJh7t/aPWx8fK5p2OdU4M7N+/jpL6NtLU0ifHndDjbtOsRBZ0cKI8roXW5MwRh1RTHeoGII20oAtLXkEIEjx/KR5zZLWCTF0NeOWv385In9oeOFcgQ5ejz6uX5V7ty8L1YeiDZ6lxtLdjTqirTeoGIJWxkBTBg3hp4Y5QJErmzgdcUVNSfu3KwIM25nhSkYo65YvngWLbmhqWjluGHiVkZJyksg0kUcpbh8klY/5UTwjMBZGaXDMAVj1BWd8zr42iVn0NHWUtYbJm5lFKbUgihERgTHbd0EWHBqe+xrl4uOthb2rrqQx1csqqhb3WwwRt2RxhtUDFHGXX9l5L/X6g27U3tvfM9Q3AZIgSf3HmHJ26exadehVMbiUqjklqgQW8EYDY1vIyk04ra15IasjDrndfD4ikV0pLABBaOBk8gPKA9tfzH2tUthpCu8ckUJ2wrGaGjijLthN2VSRHCUxyiOnt585GuXSjDYrlhGEmtUiCkYo6Ep1u0d3C6FBczFeYzimLHiYdpbc1x6dgcPbX9xUOmUQpqVUFyUcDljjWyLZDQ0xbq9i70xi+HIsTxrntrPSBxLAon2lqSEznLGGpmCMRqKQtvCwtmTyDUPvaNzzRJ6k4bdmNes2cbcG37Buq7uskTH5vs1NqgvCSV5G5NUG6ecsUamYIyGIUxBrHlyP/0DBVuaiB1O1Aqlpzc/GN5fbdJsj5JWKOWMNTIbjNEwhCmIfKFycWOfvWc7MHQ1ELdC6c3301S5mLlIFs6exLmrNsYmVDaJhNqJ/BVKnJ2pWKxcg9EwRJViiKIl18ylZ3cMJg5G3Zi1REuueZiH69KzO4YkVEadN5KARSvXYDQ8UaUYoujN93PX5n2DSqnWlUuTEGpbCUuoBC9NYUA10yJUpmCMhiEsziTXJCCecTWMsNEmgZCdVVVpbpLhtiRHlGIcUGXvqguzFMuMvEZjETTEtrXkWH35HFZfNqeopMMBhfbWHG0tuSxELJrWXBM3XT4n0sAb9dkqUbLBFIzREISlBLzWNwB4Rs2bPjRnmOckTuUcOZbntb4BWmvAc9Q+YRyd8zoivT9XnDMt9LN19/RmXiy8+lfHMCpAUuxHWJb2lQumx2Y69+b7OZYfyFLsVPh2pahM8xs7zxgcB0+5aODcLLsmmA3GaAiSYj+i6uLGlaesFfx6NH6WeZix1h8PKxaeZY1eUzBGQxDlQZrS1hKa3Lf8vu2g4XEytYZfjyaNgsiq5GgUpmCMUcO6rm5ueHDnoJ2lrSXH9RedPmifiMqCDg3Ai/Aq1Sq+PSUpMC5O0WZBpjYYEWkTkftEZJeIPCsi7xCRiSLyiIg85363u7kiIt8SkT0i8rSInBV4nWVu/nMisiz6HY1GZV1XN8vv2z7EiNvTm2f5vdsHtw9RlfAqVWE/S3yjbVI3yqxKjkaR9Qrmm8DPVfUyERkLtAJfBB5V1VUisgJYAXwBuACvLexpwDnArcA5IjIRuA6Yj7ca3Coi61X1SMayG3XE6g27Q1cd+QEd3D6E2SfWdXXXRYRuHEGjrU+UXaWcaQBpyEzBiMiJwH8EPgygqseB4yJyMfBuN+0O4DE8BXMx8EPXj3qzW/1MdnMfUdXD7nUfAd4H/CQr2Y3syKqNadwqJOq5uPotuWapyjapo62FhbMnDaYnJEUfd8Q8H1fTplJ1ebNcwZwKHAJ+ICJzgK3A1cCbVPVFAFV9UUTe6OZ3APsD5x9wY1HjQxCRq4CrAKZPn17eT2KURKEyWTh7UmiTMSi+UlohcTdilH0hKju6SWDC2DEjKvpUCrcsmRt6HZLaxEY9X6neR3FkaYMZA5wF3Kqq84CjeNuhKMLimjRmfOiA6m2qOl9V50+aNKkUeY0yElYa4a7N+2JjUUbC8sWzaA5JZ24iugBT1Df8gJKJcokL3GtvzUUq2SS7SaXtKsWQpYI5ABxQ1Sfc8X14Cud3buuD+/1SYP60wPlTgYMx40YNE7Y6yKKNqV9A6po120JzcSSmhkKlv+GvXDA9NL2gJdfMdR88PfK8pFYtWbVyKQeZbZFU9bcisl9EZqnqbuA84Bn3swxY5X4/4E5ZD3xaRO7GM/K+4rZQG4C/971NwPnAyqzkNspDMUqj2Bvd33p19/SGGjiD9AeMvIXMOLm47OqR0N6a48ZOL6q2FDtUkt2kknaVYsjai/RfgbucB+k3wEfwVk33iMjHgH3A5W7uz4D3A3uAY24uqnpYRL4KPOXmfcU3+BojIyuDK6QvjeAv5dPKUhgUl8YMG9az6Pr1OytqYwmuUAo9OcF0hdFGpgpGVbfhuZcLOS9krgKfinid7wPfL690jU0xrSmCK4Zm59LtSFBIaVpwCHDp2d75aWUppbB2Yc+i5fdur2iErm9fiVp5ldPYXWtYRbsGJckz4VOoiIL4N0mUsknbI6g5Ig6lUBYovipdk8BJLTl6juWZ0tbCS6/2Usn8RL9SHJCocMM+b5arzHJiFe2MIaTNSYlbMSR9A/uPk26sKAVUGP6+cPakooPiBpTB6N5K2Vt82ltzXPfB0weTDJNWXmFbuXI1QKsWqRWMiExQ1aNZCmOUhzTfemlzUtIaa5MiR4vtdujjy9jd08udm/cVfX41CFvRpbmOhde+nA3QqkWim1pE3ikizwDPuuM5IvKdzCUzSiKpqZZP2tiJYjw8cZGjYQWdRhstuWZuWTKXx1csClXoSecWXvtKZz5nQZo4mJuBxcDvAVR1O14KgFGDJBVW8kkbOxGmiKIovImCTc5Wb9jNpWd3xPbt8Us7FlO+slZIij0Ju47+p4w6t5wN0KpFqi2Squ6XoX/0kXfnNjKhmG+9NLETQZdq0ItUGH9S+A0cZj+4f2t3pMHTb6+x5sn9dVGDJYj/2ZPiVKC4JMO4EhP1QhoFs19E3gmoi2f5DG67ZNQeWdT7iMpCjrtZolZS16/fyYRxY+jN9w9zeV+/fmfdKRdIbxcpNhiu0pnPWZBGwXwCr+xCB17Y/i+IiFcxqk+lvvWSbpYoj01Pb34wwK1flVyTsHD2JFZv2F3x5MJykpVdpFYjdNOSqGBU9WXgygrIYpSBWvnWi4ptKSQ/oHXjHYqjnuwilSRRwYjIHcDVqtrjjtuBm1T1o1kLV+9UK0iqFr716rmAUyksnG0Z/GGk8SKd6SsXAFdJbl52Io0O0rqLRytx3qLRyKZdh6otQk2SRsE0BTKZcSUsLQI4gbTu4nIRdAln3UwrDcW4t0cD9RSbUknSKIqbgH8Tkfvc8eXAf89OpNFBJYOkai2k3N8aFnqKFs6exJqn9tddxf40mA0mnDRG3h+KyBZgEV5s0CWq+kzmktU5lWwPUa2Q8qhmZUFl1686LE7kJ0/srykbTVtLDhGGdCQohnqLTakkkdnUInKiqr7qtkTDqOWaLLWQTR2Whexn1pb7po/KMBZg76oLy/pePlGfT9DQdqpRAXq1gN8/KSkpMzh/wrgxdRubkgWlZFP/GPgAXrHu4P+E/z9yalklHGVU0l1c6WZaEL1qisJfsdSacgEvNqeYiGW/mZuRTKSCUdUPiJcf8JeqWv+BClWgUu7iaoSUj0ajZikRy0Y8sTYYVVUR+SlwdoXkMUpgJKulpBso6vm0JTHrgfbW4YW4fWohpqieSeNF2iwib1PVp5KnGtWilBshyfvkt2P1vT6DTeFJVxKzXrjwzMlDjm3VUj7SxMEsBP6fiPza9YzeISJPZy2YkT1JsTo3PLhzmEs536/c8ODOIeUeoHZLLLTkmoaUpDj3zROH9Se6f2v3YNxQowdIlps0K5gLMpfCqApJsTpRblt/PG1JzGrytUvOHLL6OHfVxtg+zqOhilwtkbiCUdUXgJPxekdfBJzsxow6pxwFjUqp8l8p2lqGd0tMUqqjoYpcLZGmZObf4TWpPxk4Ba/X9JezFszInqSymWFdCAvHa/nGu/6i4d0Sk5TqaKgiV0ukscFcAbxNVa9T1euABVj5hlFBUtnM6y86nVxB69VckwzeuOu6ummqUdsLhKdJJCnVqKxoy5YujTQ2mOeB8cCf3PE44Ndp30BEmoEtQLeLrZkJ3A1MBH4J/LWqHheRccAP8VzivweWqOrz7jVWAh/DK9X5GVXdkPb9jXiSvE9vGD9m0ObiR7z6HqaVa3fUVMh/kPbWXKw3KGo8KivasqVLI42CeQ3YKSKP4AU1vhf4VxH5FoCqfibh/KvxSmye6I6/DtysqneLyHfxFMet7vcRVf1zEVnq5i0RkbcAS4HTgSnA/xaRv1DV2tz41wFpYl9Wrn2a3oKQ/9f6Xj+uZdtLrlm48MzJsS74KKVqNpjykmaL9FPgi8Am4DHgS8A/46UQbI07UUSmAhcC33PHgpc06Wdm3wF0uscXu2Pc8+e5+RcDd6vqa6q6F6939dtTyG2EkOSG/fK6HVyzZtsw5QKeN+WaNdsiu0JWk2CF/tWXzWHTrkMllcswG0x5SZNNfUfSnBhuAT4PnOCOTwZ6VLXPHR/Aq/WL+73fvWefiLzi5ncAmwOvGTxnEBG5CrgKYPr06SMQeXQS7ItcSPDGuytF+cpaUy7wegtbv/XqtWu2hc5LWomMhkr+tUSaFUxJiMgHgJdUNbjKCbMIasJzcee8PqB6m6rOV9X5kyaZQS5IcNUSxcGeXlZv2F2TyYg+QnxYf1B5lLoSSdsvykhHlpXpzgUuEpH34xmJT8Rb0bSJyBi3ipkKHHTzDwDTgAMiMgY4CTgcGPcJnmOkII29ZEpbS83bGRTo+rvzI7doQeUxkpWI5R+Vj8xWMKq6UlWnquoMPCPtRlW9Es+Wc5mbtgx4wD1e745xz29Ur1jNemCpiIxzHqjTgCezkns0kqQ4/Buv1u0MflpCmra35VqJ1Fop0nojcgUjIg8SU75DVS8q8T2/ANwtIjcCXcDtbvx24Ecisgdv5bLUvc9OEbkHeAboAz5lHqThxHmGkjKfx+e875laTmDMNcugAkmbPT7SlUitlSKtR+Iq2v2le3gJ8B+AO93xFcDzqvrF7MUrjVqoaFdJwqrL+YWS/Fq492/tjlUcfrU9eP3GbWvN8cqxPMP9SZXnliVzK35TR23FgsZkw6Poinaq+n/ciV9V1WCz+wdF5P9mIGNFqYeU/LQyhtlY/K+N7p5e7ty8j5ZcE+2tucgERt+TFNxm9BzL14TRt1mkKn8bi4kZOWmMvJNE5FRV/Q2As4PUtZumUkvfkSixYmRM8w/fmx+gr1/JNUtkVX//PWpti9Svyrqu7oormWqUIh1tpDHyXgs8JiKPichjeEbaazKVKmMq0bNopHVFbnhwZ2oZ0/7D5wc0tmVIs0hVlIsI/NWC6bE1ZapRkyWNMdmIJ025hp/jeW6udj+z6j0XqBJL35EosXVd3ZFbme6e3mFejXI1OatWXpGqV/TpinOmRX6OLJvWRWExMSMnTbmGVmA58GlV3Q5Md0F0dUslwsFHosTibiSBYasiYEh1uXqkN9/Ppl2HBg3NYVTD9tE5r4PHVyxi76oLeXzFIlMuRZJmi/QD4DjwDnd8ALgxM4kqQCWWviNRYnE3Ulg1ts/esz0yND5IrknINceXV0h6PksO9vTSOa8jUlGa7aP+SKNg3qyq/wDkAVS1l/Dw/bqhEkvfkSixYm+kftXBFU3UH6ZZhNWXz2HJ26ZFzPBsIVm0de1oa+GvFkxnwtj4bdyUIgLpjPogjRfpuIi04L48ReTNeCUc6pqsw8FH0kpk+eJZQ6r5g7ey6BtQkswkfvJWYbMwf+tx/9ZoQ2kWJphmkVQBfEEFUsmmdUa2pFEw1wM/B6aJyF14OUYfyVKo0cKIGnkV3uyaXgEEpwWLRM294RcV9xL1q0bmQjWLMKAaeh0sH2h0kKZcwy9EZCteqUwBrlbVlzOXbBSSNrZl9Ybd5AcK2oUM6GA702Lwi0St6+qmp7e05u4joSMmiXJANbPe2UZtkMaL9Kiq/l5VH1bVh1T1ZRF5tBLCjTbSuq6jbsh+1aLd0f7rV9rFC69ve9oiSixEjRujh7hkx/FAK3CKiLTzumH3RLzSlUaANFufKMXhx7b4557UkgtdbfiBcMWuZKpRIKq9Ncd1H/S2Ztev3xk6p0bL+RplJG6L9F/wInan4JXG9BXMq8C3M5arrki79YkKPfdjW/xzc81CrkmGbZN8pVKpgLjn3fallBKZrWPHDH72VyK2ZlHjxughcoukqt9U1ZnA51T1VFWd6X7mqOr/rKCMNU/arU+Y+7XQ4wOeq7hQuVQDPzQ/KlI4qm8SlKe6nFH/pImDGRCRNv9ARNpF5JMZylR3pI3aDYu/qb4aicbP/wmT+5Ylc9l23fmpguIsrqVxSeOm/riqDm6JVPWIiHwc+E52YtUXxWTdFrpfa7FCv0+wJ3OU2zhNaUqLa2lc0iiYJhERV77Sb6Q2Nlux6ouwm8y3q5y7amPkzbSuq5ujr/UNG68lkvJ/KlVdzqhP0iiYDcA9rkmaAp/AC7wzHMGbzF+NBAs+Lb93+5B5EF6FrhZJYycx5WFEkcYG8wVgI/C3wKeAR/F6HRkB/KzbMMNnfkCHuWor3RnxTSeMTUwgK3ze7CTGSEkTyTuA19r11uzFqX+iomV7evNDYmVKNe6GeZ3S8PIf85G2Ir/GbD2UETXqi7hAu3tU9UMisoPwRmdnZirZKKQcW6IxTRDS1TWRftVEg6xtdYxyE7eCudr9ruviUpUmqrB2k1CWLVEpyiXI1y45w1YpRsWIbFtSz1Szbcm6ru7QUgtZ1FkpFmu3YWRFVNuSSCOviPxBRF6N+knxhtNEZJOIPCsiO0Xkajc+UUQeEZHn3O92Ny4i8i0R2SMiT4vIWYHXWubmPyciy6LesxbonNfB6svmDAlK848rQVzh7DTxNtbJ0CgncX2RTgAQka8AvwV+hGdjvBI4IcVr9wGfVdVfisgJwFYReQT4MPCoqq4SkRXACjxP1QV4xcVPA87BMyqfIyITgeuA+Xi2oK0isl5Vj5TweStClC0jqjlaORmIWZHGKR+wToZG+Unjpl6sqt9R1T+o6quqeitwadJJqvqiqv7SPf4D8CzQAVwM3OGm3QF0uscXAz9Uj81Am4hMBhYDj6jqYadUHgHeV8RnrAnCwu2vXDC9LN0AfDraWmLjVuKSJNd1dfPZe7Zn3s7FaCzSBNr1i8iVwN14X7hXAEVZK0VkBjAPeAJ4k6q+CJ4SEpE3umkdwP7AaQfcWNR43RG2spn/ZxOHBOiViu8N2vLCYe7cvC90TkdbS6grGrzVVZQCsk6GRqmkUTD/Cfim+1HgcTeWChF5A3A/cI2qvirRy/SwJzRmvPB9rgKuApg+fXpa8TKhmHgSf3z5vdtLzqBuzTXx965oedRqQ4CFsyeFboHG55piPVyW9WyUSppAu+fxti9FIyI5POVyl6qudcO/E5HJbvUyGXjJjR8AgiXvpwIH3fi7C8YfC5HzNuA28LxIpchbDkqxY4SVyCyGY/kBrlmzLXYlpMCmXYdCt0Bpi3EbRrGkKZn5FyLyqIj8yh2fKSJfTnGeALcDz6rqPwaeWg/4nqBlwAOB8b9x3qQFwCtuK7UBON+ViWgHzndjNUkxHR19j025sqnjXqetJVf0VqdZxDoZGiMijZH3fwEreb0v0tPA0hTnnQv8NbBIRLa5n/cDq4D3ishzwHvdMcDPgN8Ae9x7ftK932Hgq8BT7ucrbqwmSVsbJti7uhKIRG91WnJNofVabvrQHFMuxohIY4NpVdUnC2wniTUGVPVfiW7Qdl7IfMVLpgx7re8D308WtfqkrQ0T1tw+S3qO5bnug6eH2nr6BpQlb5vKpl2HLMLXKCtpFMzLrtmaXw/mMuDFTKWqY9IUYIprbp8VU9pa6JzXwQ0P7hz23vl+ZdOuQxbla5SdNArmU3jG09ki0g3sxQu2G/WUkl2cpgBTpeNKck3CseN9zFzxcGRgn7mijSyItcGISBMwX1XfA0wCZqvqu1T1hYpIV0WCNhK/77NfozbpPF+5tLXmOPpaH9eu2TYk7H6kN3NHWwvPr7qQ51ddyC1L5sZG6La15EDgyLF8bNSwuaKNLIhVMK4WzKfd46MuIrchKMYb5FOolI4cy9PTmx9UUNeu2caMFQ/TlBCyn0RQQXXO64hMDxBgwrgxiYmW5oo2siKNF+kREfmcS16c6P9kLlmVSesNCpJUpc6/zUfa16hJZMhKKq4tSNJqqaOtxVzRRmakscF81P0OengUOLX84tQOabxBhTaaSrmc+1VZuXYHW144zKZdh+ju6R2WOOmvSuKC7wTMsGtkSppI3pmVEKTWSPIGhUXsZpEdHUVvvn9IzpGfU6F4q5KgYfnaNdtC5TK7i5E1iQrG9aj+JPAuvP/ffwG+q6p/yli2qpLkDQrbDpVTuTQJFJs94CuX4Kqkc14HW144zF2b94WucAwjS9JskX4I/AH4H+74CrzaMJdnJVStEFejNmu3bqlmmjC5buw8YzBr2wLpjEqSRsHMUtU5geNNIrI9K4HqhaxtLie15CI7FMQRte2xgt5GNUjjRepyyYcAiMg5eCUbGprli2cl9hkaCUeP95FrKu4dbNtj1BppVjDn4GU5+xbF6cCzfjuTRmpfUug1ytKgm+9X2ltztI4dMxi098c/9Q3JI8o1CW8YP4aeY3nb9hg1SRoFU3flKctFYVRu8AavhNeo51ierr87P1QeUyhGPZDGTT3q0wLCKHRDhyUnZu2SLrSnmB3FqDfS2GAakmJ6R08Y21x2e4zZU4zRgCmYCIpxQx873s/NS+bS3jq88X0pWCU5Y7RgCiaCYqJcFa+A1B//lFiHawjtrTmrJGeMakzBRLB88ayiehYdOZYPLdwdtXVqyTVz3QdPH9YryVYuxmgijRepIQlLFVg4e9KwkPskgnOjcoVMoRijFVMwMUR5bQqVTFp3dViukGGMZmyLVCQ3dp7BzUvmDjazLzYWxkpTGo2EKZgS6JzXweMrFtERE80bVcbSSiQYjYQpmBEQtRoR4KYPzQn1EFlsi9FImIIZAXGlKjvndZiHyGh46sbIKyLvA74JNAPfU9VVCadkTlLVOwvtNxqdulAwItIMfBuv1ewB4CkRWa+qz1RTrjQ9kAyjkakLBQO8Hdijqr8BEJG7gYuBkhXMwMAAAE1NI9sl2irFMKKpFxtMB7A/cHzAjQ0iIleJyBYR2XLo0KHEF2xqauL48ePlldIwjCHUi4IJ8/kO8RCr6m2qOl9V50+aNCnVi44fP74cshmGEUG9KJgDwLTA8VTgYJVkMQwjJfWiYJ4CThORmSIyFlgKrK+yTIZhJFAXRl5V7RORTwMb8NzU31fVnVUWyzCMBOpCwQCo6s+An1VbDsMw0lMvWyTDMOoQUzCGYWSGKRjDMDLDFIxhGJlhCsYwjMwwBWMYRmaYgjEMIzNMwRiGkRmmYAzDyAxTMIZhZIYpGMMwMqNucpGqybqubiuLaRglYAomgXVd3UMKe3f39LJy7Q7AWr4aRhK2RUpg9YbdQ7oGAPTm+1m9YXeVJDKM+sEUTAJRzdWsBaxhJGMKJoG45mqGYcRjCiaB5YtnWQtYwygRM/ImYM3VDKN0TMGkwJqrGUZp2BbJMIzMMAVjGEZmmIIxDCMzRFWTZ9UZInIIeCHF1FOAlzMWJwtM7spicifzZ6o6rGfzqPLwJvwAAAW4SURBVFQwaRGRLao6v9pyFIvJXVlM7tKxLZJhGJlhCsYwjMxodAVzW7UFKBGTu7KY3CXS0DYYwzCypdFXMIZhZIgpGMMwMqNhFYyIvE9EdovIHhFZUWVZponIJhF5VkR2isjVbnyiiDwiIs+53+1uXETkW072p0XkrMBrLXPznxORZRWSv1lEukTkIXc8U0SecDKsEZGxbnycO97jnp8ReI2Vbny3iCyugMxtInKfiOxy1/0d9XC9ReRa9z/yKxH5iYiMr+nrraoN9wM0A78GTgXGAtuBt1RRnsnAWe7xCcC/A28B/gFY4cZXAF93j98P/DMgwALgCTc+EfiN+93uHrdXQP7/BvwYeMgd3wMsdY+/C/yte/xJ4Lvu8VJgjXv8Fvc3GAfMdH+b5oxlvgP4z+7xWKCt1q830AHsBVoC1/nDtXy9q3JDVfsHeAewIXC8ElhZbbkC8jwAvBfYDUx2Y5OB3e7xPwFXBObvds9fAfxTYHzIvIxknQo8CiwCHnI34cvAmMJrDWwA3uEej3HzpPD6B+dlJPOJ7kaVgvGavt5Owex3Cm2Mu96La/l6N+oWyf9D+RxwY1XHLWPnAU8Ab1LVFwHc7ze6aVHyV+Nz3QJ8HhhwxycDParaFyLDoHzu+Vfc/ErLfSpwCPiB29p9T0QmUOPXW1W7gW8A+4AX8a7fVmr4ejeqgpGQsar760XkDcD9wDWq+mrc1JAxjRnPBBH5APCSqm4NDsfIUBNy432bnwXcqqrzgKN4W6IoakJuZxO6GG9bMwWYAFwQI0PV5W5UBXMAmBY4ngocrJIsAIhIDk+53KWqa93w70Rksnt+MvCSG4+Sv9Kf61zgIhF5Hrgbb5t0C9AmIn4xs6AMg/K5508CDldB7gPAAVV9wh3fh6dwav16vwfYq6qHVDUPrAXeSQ1f70ZVME8Bpznr+1g8A9j6agkjIgLcDjyrqv8YeGo94HsmluHZZvzxv3HejQXAK25JvwE4X0Ta3bfd+W4sE1R1papOVdUZeNdwo6peCWwCLouQ2/88l7n56saXOq/HTOA04MkM5f4tsF9E/MLK5wHPUOPXG29rtEBEWt3/jC937V7vrAxStf6D5xn4dzwL+peqLMu78JaoTwPb3M/78fbLjwLPud8T3XwBvu1k3wHMD7zWR4E97ucjFfwM7+Z1L9Kp7h92D3AvMM6Nj3fHe9zzpwbO/5L7PLuBCyog71xgi7vm6/C8QDV/vYEbgF3Ar4Af4XmCavZ6W6qAYRiZ0ahbJMMwKoApGMMwMsMUjGEYmWEKxjCMzDAFYxhGZpiCMVLjMpA/WW05AETk36otg5GMuamN1Lg8qYdU9a0hzzWran8FZKjI+xjlwVYwRjGsAt4sIttEZLWIvFu8OjY/BnaIyAwR+ZU/WUQ+JyLXu8dvFpGfi8hWEfkXEZld+OIicr2I/EhENrraJh9340Pex439MXDe50Vkh4hsF5FVad/PyJ4xyVMMY5AVwFtVdS54Nz7wdje2N1jQKITbgE+o6nMicg7wHbzcpULOxKu5MgHoEpGH3fjg+wQni8gFQCdwjqoeE5GJRb6fkSGmYIyR8mThTV+IyxJ/J3Cvl0IDeCHuYTygqr1Ar4hswlMsPTHv8x7gB6p6DEBVDxf5fkaGmIIxRsrRwOM+hm67x7vfTXg1S+ameL1Co6B/fLRwokNCzinm/YwMMRuMUQx/wCvpGcXvgDeKyMkiMg74AIB6tW32isjlMFjjdk7Ea1zs6syejJdA+VSCTL8APioire61Jxb5fkaGmIIxUqOqvwceF6/g9OqQ5/PAV/Cq8T2El/XrcyXwMRHZDuzEK5wUxpPAw8Bm4KuqGlunRFV/jld+YIuIbAM+V+T7GRlibmqjZnAepz+q6jeqLYtRHmwFYxhGZtgKxjCMzLAVjGEYmWEKxjCMzDAFYxhGZpiCMQwjM0zBGIaRGf8fw/HOw1xv710AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#最小二乘线性回归\n",
    "# 1. 线性回归\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_train, y_train)\n",
    "y_test_pred_lr = lr.predict(X_test)\n",
    "y_train_pred_lr = lr.predict(X_train)\n",
    "# 权重系数\n",
    "fs = pd.DataFrame({\"columns\":list(feat_names), \"coef\": list((lr.coef_.T))})\n",
    "print(fs.sort_values(by='coef', ascending=False))\n",
    "\n",
    "# RMSE 评估模型性能\n",
    "print('The RMSE of LinearRegression on test is %.6f'  % RMSE(y_test, y_test_pred_lr))\n",
    "print('The RMSE of LinearRegression on train is %.6f'  % RMSE(y_train, y_train_pred_lr))\n",
    "\n",
    "# 观察残差分布，看是否符合：噪声为0均值的高斯分布\n",
    "f,ax = plt.subplots(figsize=(7,5))\n",
    "f.tight_layout()\n",
    "ax.hist(y_train-y_train_pred_lr, bins=40, label='Residuals Linear', color='b', alpha=0.5)\n",
    "ax.set_title(\"Histogram of Resicuals\")\n",
    "ax.legend(loc='best')\n",
    "plt.show()\n",
    "# 观察预测值与真值的散点图\n",
    "# plt.clf\n",
    "plt.figure(figsize=(4,3))\n",
    "plt.scatter(y_train, y_train_pred_lr)\n",
    "plt.plot([0,1],[0,1], '--k')\n",
    "plt.axis('tight')\n",
    "plt.xlabel('true price')\n",
    "plt.ylabel('predicted price')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最小二乘线性回归模型拟合的结果显示：\n",
    "- 残差基本上符合均值为0的高斯正态分布，说明本数据集训练出的模型比较合适\n",
    "- 从散点图上看出，真实值到训练值基本上在y=x直线附近，离散较远的样本在较小到较大的值之间都有分布\n",
    "- 系数的值非常大，因为特征强相关，OLS性能并不好\n",
    "- 适合特征数量不太多，而且相关性不强的时候，适合使用OLS模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The RMSE of RidgeCV on test is 796.615615\n",
      "The RMSE of RidgeCV on train is 746.549433\n",
      "[0. 0. 0.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEGCAYAAABVSfMhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5RV5Znn8e9TV67FtbhfjSiCxgslYEyMLYmikw7GSNSI0Bm7nU4nM53u6Z6Y6bXaNcmsmaTT0+kx3XHG7jiRwqBEYzSJkUZNYl8ooFBRwAuoHKqggIICqqCgrs/8sd8qD+WpK3XOPlX1+6x11tnn3e/e73O2WE/tvd96trk7IiIimZITdwAiIjK0KPGIiEhGKfGIiEhGKfGIiEhGKfGIiEhG5cUdQLabOHGiz5kzJ+4wREQGlO3btx919+JU65R4ujFnzhzKy8vjDkNEZEAxs0Rn63SpTUREMkqJR0REMkqJR0REMkqJR0REMkqJR0REMkqJR0REMkqJR0REMkqJR0REPuRvX3iHLe8dS8u+lXhEROQc71Wf4m9f2MPW92vSsn8lHhEROce6sv3k5Rh3LJ6Zlv0r8YiISLv6xmZ+sr2C5ZdOYdLoYWkZQ4lHRETa/XzHQerONnPP0tlpG0OJR0REAHB31m5OcPHk0SyeOz5t4yjxiIgIAK9WnGDXwVpWXTMbM0vbOEo8IiICwLrNCUYV5vG5K6endRwlHhER4dipBn7xehW3XTWdUYXpfVSbEo+IiLChvJLGllZWpXFSQRslHhGRIa6l1XlsS4Ilc8dz0eTRaR9PiUdEZIj7zdtHqDx+htXXzMnIeEo8IiJDXGlZguLRhdy4cHJGxlPiEREZwhLHTvPbd6q5a/Es8nMzkxKUeEREhrAfb9lPjhlfXDwrY2Mq8YiIDFFnm1p4oryCGxdMZsqY9NRlS0WJR0RkiPrF61WcqG9Ka122VNKWeMzsETM7YmY7k9pWmtkuM2s1s5IO/b9hZnvN7G0zuympfXlo22tm9ye1zzWzLWa2x8yeMLOC0F4YPu8N6+d0N4aIyFBUWpbgI8UjueYjEzI6bjrPeH4ELO/QthO4DXg5udHMFgB3AgvDNj8ws1wzywX+HrgZWADcFfoCfAf4nrvPA44D94b2e4Hj7n4h8L3Qr9Mx+u3biogMIK9XnmBHxQnuWZreumyppC3xuPvLQE2Htjfd/e0U3VcAj7t7g7u/D+wFFofXXnd/z90bgceBFRYdpRuAJ8P2jwK3Ju3r0bD8JLAs9O9sDBGRIad0c4IRBbnctmhGxsfOlns804GKpM+Voa2z9gnACXdv7tB+zr7C+pOhf2f7+hAzu8/Mys2svLq6+jy+lohI9jlR38izOw5y65XTKRqWn/HxsyXxpDrP8z6092VfH250f9jdS9y9pLi4OFUXEZEB6yfllTQ0t7JqSWYnFbTJlsRTCSQ/3HsGcLCL9qPAWDPL69B+zr7C+jFEl/w625eIyJDR2uqs25KgZPY4FkwriiWGbEk8zwJ3hhlpc4F5wFZgGzAvzGArIJoc8Ky7O/Br4Paw/RrgmaR9rQnLtwMvhf6djSEiMmT8896jJI7Vc8818ZztAKTtoQtmth64HphoZpXAA0RnHt8HioFfmtlr7n6Tu+8ysw3AbqAZ+Iq7t4T9fBXYCOQCj7j7rjDE14HHzey/A68CPwztPwRKzWxvGO9OgK7GEBEZKko372PiqAKWXzolthgsOhmQzpSUlHh5eXncYYiInLeKmnqu++6v+aPrP8Kf3zQ/rWOZ2XZ3L0m1LlsutYmISJqt37ofA74Y06SCNko8IiJDQENzC09sq2DZJZOZPnZ4rLEo8YiIDAG/euMQx043ZrwuWypKPCIiQ0BpWYK5E0fy8Qsnxh2KEo+IyGC36+BJtieOc/eSWeTkZLYuWypKPCIig9y6sgTD8nNYuWhm950zQIlHRGQQO3mmiZ+9epAVl09nzIjM12VLRYlHRGQQe2p7JWeaWmKtVNCREo+IyCDl7qwrS3DFzLFcOn1M3OG0U+IRERmk/u3dY7x39DSrs+hsB5R4REQGrbWb9zFuRD63XDY17lDOocQjIjIIVZ08w6bdh/nC1TMZlp8bdzjnUOIRERmE1m/Zj0NsD3vrihKPiMgg09jcyvptFfzOxZOYOX5E3OF8iBKPiMggs3HXIarrGrKiLlsqSjwiIoNMaVmCmeOHc91FxXGHkpISj4jIIPL2oTq2vl/DqiWzyc2CumypKPGIiAwipWX7KMjLYWVJdtRlS0WJR0RkkKg728TTrxzgMx+dyviRBXGH0yklHhGRQeJnrx7gdGMLq6+ZE3coXVLiEREZBNydtZsTXDZ9DJfPyJ66bKko8YiIDAJb3q9hz5FT3LN0NmbZOamgjRKPiMggUFqWYMzwfH738mlxh9ItJR4RkQHuSO1ZNu48xMpFMxhekF112VJR4hERGeDWb62gudW5O0srFXSkxCMiMoA1tbTy460JrruomLkTR8YdTo8o8YiIDGAv7D7M4drsrcuWihKPiMgAVlqWYPrY4dwwf1LcofRY2hKPmT1iZkfMbGdS23gz22Rme8L7uNBuZvagme01s9fN7KqkbdaE/nvMbE1S+yIzeyNs86CF+YN9GUNEZCDae6SOf3v3GF9cMitr67Klks4znh8Byzu03Q+86O7zgBfDZ4CbgXnhdR/wEERJBHgAWAIsBh5oSyShz31J2y3vyxgiIgPVurL95Ocad1ydvXXZUklb4nH3l4GaDs0rgEfD8qPArUntaz1SBow1s6nATcAmd69x9+PAJmB5WFfk7pvd3YG1HfbVmzFERAac0w3NPLW9klsum8rEUYVxh9Mrmb7HM9ndqwDCe9tFyelARVK/ytDWVXtliva+jPEhZnafmZWbWXl1dXWvvqCISCY889pB6hqaWX3NwJlU0CZbJhekujjpfWjvyxgfbnR/2N1L3L2kuDg7H6QkIkNXVJdtH5dMLeKqWeO67Z9tMp14Drdd3grvR0J7JZB8kXIGcLCb9hkp2vsyhojIgLI9cZy3DtUNiLpsqWQ68TwLtM1MWwM8k9S+Osw8WwqcDJfJNgI3mtm4MKngRmBjWFdnZkvDbLbVHfbVmzFERAaU0rIEowvzuPXK7K/LlkpeunZsZuuB64GJZlZJNDvt28AGM7sX2A+sDN2fA24B9gL1wJcA3L3GzL4FbAv9vunubRMWvkw0c2448KvwordjiIgMJEdPNfDcG1XcvWQ2IwrS9iM8rdIWtbvf1cmqZSn6OvCVTvbzCPBIivZy4NIU7cd6O4aIyEDxxLYKmlqcVQOoUkFH2TK5QEREutHS6jxWluDaCydw4aRRcYfTZ0o8IiIDxEtvHeHgybMDqi5bKko8IiIDxNrN+5hSNIxPXTI57lDOixKPiMgA8P7R0/zznqN8ccks8nIH9o/ugR29iMgQ8VhZgrwc484BVpctFSUeEZEsd6axhQ3lFdx06RQmFQ2LO5zzpsQjIpLlfr7jILVnmwf8pII2SjwiIlnM3Vlbto+LJo9iydzxcYfTL5R4RESy2GsVJ9h5oHbA1mVLRYlHRCSLlZYlGFmQy61XpnyKy4CkxCMikqVqTjfyi9eruO2qGYwelh93OP1GiUdEJEttKK+gsbl1QNdlS0WJR0QkC7W0Oo9tSbB47ngunjI67nD6lRKPiEgWevmdaipqzgzIR1t3R4lHRCQLrd28j+LRhdy4YErcofQ7JR4RkSyz/1g9v3mnmruunklB3uD7MT34vpGIyAD32NYEOWbctWRW3KGkhRKPiEgWOdvUwoZtFXz6kslMHTM87nDSQolHRCSL/PL1Ko7XN3HPIJxU0EaJR0Qki5SWJbigeCQf+8iEuENJGyUeEZEs8UblSV6rODGo6rKlosQjIpIlSsv2MTw/l9uumhF3KGmlxCMikgVO1jfxzGsHufXK6YwZPnjqsqWixCMikgV+sr2ChuZWVi0dnFOokynxiIjErLXVWVeWYNHscSycNibucNJOiUdEJGb/svco+47VD8q6bKko8YiIxKy0LMGEkQUsv3Tw1WVLRYlHRCRGB06c4cU3D3PH1TMpzMuNO5yMUOIREYnRj7ckAPjiIK3LlkqPE4+ZfdzMvhSWi81sbl8HNbM/NrOdZrbLzL4W2sab2SYz2xPex4V2M7MHzWyvmb1uZlcl7WdN6L/HzNYktS8yszfCNg9a+EuszsYQEYlDQ3MLT2yr4Ib5k5kxbkTc4WRMjxKPmT0AfB34RmjKB9b1ZUAzuxT4A2AxcDnwGTObB9wPvOju84AXw2eAm4F54XUf8FDYz3jgAWBJ2NcDSYnkodC3bbvlob2zMUREMu75nYc4eqpxUNdlS6WnZzyfAz4LnAZw94NAX5/FeglQ5u717t4M/DbsfwXwaOjzKHBrWF4BrPVIGTDWzKYCNwGb3L3G3Y8Dm4DlYV2Ru292dwfWdthXqjFERDKudHOCORNG8IkLJ8YdSkb1NPE0hh/iDmBmI89jzJ3AdWY2wcxGALcAM4HJ7l4FEN4nhf7TgYqk7StDW1ftlSna6WKMc5jZfWZWbmbl1dXVff6iIiKd2X2wlvLEcVYtnU1OzuCty5ZKTxPPBjP7v0RnG38AvAD8Q18GdPc3ge8QnaE8D+wAmrvYJNV/Ee9De29ifNjdS9y9pLi4uDebioj0SGlZgsK8HG5fNLjrsqXSo8Tj7n8NPAk8BVwM/KW7f7+vg7r7D939Kne/DqgB9gCHw2UywvuR0L2S6IyozQzgYDftM1K008UYIiIZU3u2iZ+9eoAVV0xj7IiCuMPJuJ5OLhgJvOTuf050pjPczPpcxc7MJoX3WcBtwHrgWaBtZtoa4Jmw/CywOsxuWwqcDJfJNgI3mtm4MKngRmBjWFdnZkvDbLbVHfaVagwRkYz56fZKzjS1cM/SOXGHEou8HvZ7GfhE+AH/AlAO3AHc3cdxnzKzCUAT8BV3P25m3ya6pHcvsB9YGfo+R3QfaC9QD3wJwN1rzOxbwLbQ75vuXhOWvwz8CBgO/Cq8ADobQ0QkI9yd0rIEl88cy2UzBn9dtlR6mnjM3evDD+zvu/tfmdmrfR3U3T+Rou0YsCxFuwNf6WQ/jwCPpGgvBy7t6RgiIpmy+d1jvFt9mr9eeXncocSmp5MLzMyuITrD+WVo62nSEhGRoLQswdgR+Xzmo1PjDiU2PU08f0z0x5Y/dfddoWrBS+kLS0Rk8Kk6eYZ/2n2YO0pmMix/aNRlS6WnZy31QCtwl5mtIpqy3KspyiIiQ936rRW0unP3kqFVqaCjniaex4A/I/rjz9b0hSMiMjg1tbSyfut+rr+omFkThk5dtlR6mniq3f3naY1ERGQQ27jrENV1DUOuLlsqPU08D5jZPxIV1mxoa3T3n6YlKhGRQaZ0c4IZ44bzyYtSVuoaUnqaeL4EzCeqSt12qc0BJR4RkW68c7iOLe/XcP/N88kdYnXZUulp4rnc3S9LayQiIoNU6eYEBXk5fKFkZvedh4CeTqcuM7MFaY1ERGQQOtXQzE9fqeQzl01l/MihV5ctlZ6e8XwcWGNm7xPd4zGiogIfTVtkIiKDwNOvHuB0Y4smFSTpaeJZ3n0XERFJ5u6s25zg0ulFXDFzbNzhZI0eJR53T6Q7EBGRwWbr+zW8fbiO73z+MqJi+QI9v8cjIiK9VFqWoGhYHp+9fHr3nYcQJR4RkTQ4UneW53ceYmXJTIYXDN26bKko8YiIpMHjWytobnVWLdWkgo6UeERE+llzSys/3rKfT8ybyNyJI+MOJ+so8YiI9LMX3jzCodqz3KOznZSUeERE+llp2T6mjRnGDfNVly0VJR4RkX6098gp/nXvMe5eOpu8XP2ITUVHRUSkHz22JUF+rqkuWxeUeERE+kl9YzNPbq/k5kunUjy6MO5wspYSj4hIP3nmtYPUnW1mteqydUmJR0SkH7g7pZsTzJ8ymkWzx8UdTlZT4hER6Qev7D/O7qpa7rlmtuqydUOJR0SkH5RuTjC6MI9br1Bdtu4o8YiInKejpxp47o1DfH7RDEYW9vRpM0OXEo+IyHl6YlsFjS2trFo6K+5QBgQlHhGR89DS6vx4y34+9pEJXDhpdNzhDAixJB4z+xMz22VmO81svZkNM7O5ZrbFzPaY2RNmVhD6FobPe8P6OUn7+UZof9vMbkpqXx7a9prZ/UntKccQEemrX791hAMnzqguWy9kPPGY2XTgPwEl7n4pkAvcCXwH+J67zwOOA/eGTe4Fjrv7hcD3Qj/MbEHYbiHRo7l/YGa5ZpYL/D1wM7AAuCv0pYsxRET6ZG1ZgslFhXxqweS4Qxkw4rrUlgcMN7M8YARQBdwAPBnWPwrcGpZXhM+E9cssmqu4Anjc3Rvc/X1gL7A4vPa6+3vu3gg8DqwI23Q2hohIr+07epqX36nmrsWzyFddth7L+JFy9wPAXwP7iRLOSWA7cMLdm0O3SqBtTuJ0oCJs2xz6T0hu77BNZ+0TuhjjHGZ2n5mVm1l5dXV137+siAxqj21JkJdj3LVYkwp6I45LbeOIzlbmAtOAkUSXxTrytk06Wddf7R9udH/Y3UvcvaS4uDhVFxEZ4s40trChvJKbFk5hctGwuMMZUOI4N/wU8L67V7t7E/BT4GPA2HDpDWAGcDAsVwIzAcL6MUBNcnuHbTprP9rFGCIivfLz1w9y8kyTHm3dB3Eknv3AUjMbEe67LAN2A78Gbg991gDPhOVnw2fC+pfc3UP7nWHW21xgHrAV2AbMCzPYCogmIDwbtulsDBGRXllXlmDepFEsvWB83KEMOHHc49lCdIP/FeCNEMPDwNeBPzWzvUT3Y34YNvkhMCG0/ylwf9jPLmADUdJ6HviKu7eEezhfBTYCbwIbQl+6GENEpMd2VJzg9cqTqsvWRxadCEhnSkpKvLy8PO4wRCSL/OcNO3h+ZxVl/3UZo4flxx1OVjKz7e5ekmqd5v+JiPTC8dON/Pz1g3zuqulKOn2kxCMi0gs/2V5BY3OrJhWcByUeEZEeam111pXtZ/Gc8cyfUhR3OAOWEo+ISA/9dk81+2vquUePtj4vSjwiIj20bnOCiaMKuWnhlLhDGdCUeEREeqCipp6X3j7CXYtnUpCnH53nQ0dPRKQHHtuynxwzvrhEddnOlxKPiEg3zja1sKG8gk9dMompY4bHHc6Ap8QjItKN596oouZ0I/csnRN3KIOCEo+ISDdKyxJcUDySay+cEHcog4ISj4hIF3YeOMmr+0+waonqsvUXJR4RkS6Ubk4wPD+Xzy+aEXcog4YSj4hIJ07WN/HMjgPceuU0xgxXXbb+osQjItKJJ1+p5GyT6rL1NyUeEZEUorpsCa6aNZaF08bEHc6gosQjIpLCv757lPePnmb1NXPiDmXQUeIREUmhdHOC8SMLuPky1WXrb0o8IiIdHDhxhhfePMwdV8+kMC837nAGHSUeEZEO1m/ZjwN3qy5bWijxiIgkaWxu5fFt+1k2fxIzxo2IO5xBSYlHRCTJr3ZWcfRUo6ZQp5ESj4hIknVlCWZPGMF184rjDmXQUuIREQnerKpl277jrFoym5wc1WVLFyUeEZFgXVmCwrwcblddtrRS4hERAWrPNvH0qwf43cunMW5kQdzhDGpKPCIiwNOvHKC+sYXV12hSQbop8YjIkOfulJYluHzGGD46Y2zc4Qx6SjwiMuRtfu8Ye4+c0hTqDMl44jGzi83staRXrZl9zczGm9kmM9sT3seF/mZmD5rZXjN73cyuStrXmtB/j5mtSWpfZGZvhG0etPDYwM7GEJGhbV1ZgrEj8vndy6fFHcqQkPHE4+5vu/sV7n4FsAioB54G7gdedPd5wIvhM8DNwLzwug94CKIkAjwALAEWAw8kJZKHQt+27ZaH9s7GEJEh6nDtWTbuOswXSmYyLF912TIh7ktty4B33T0BrAAeDe2PAreG5RXAWo+UAWPNbCpwE7DJ3Wvc/TiwCVge1hW5+2Z3d2Bth32lGkNEhqgfb9lPq7vqsmVQ3InnTmB9WJ7s7lUA4X1SaJ8OVCRtUxnaumqvTNHe1RgiMgQ1tbSyfut+PnlRMbMnjIw7nCEjtsRjZgXAZ4GfdNc1RZv3ob03sd1nZuVmVl5dXd2bTUVkANm0+zBH6hq4R5MKMirOM56bgVfc/XD4fDhcJiO8HwntlcDMpO1mAAe7aZ+Ror2rMc7h7g+7e4m7lxQXq16TyGC1dvM+po8dzvUX6+JHJsWZeO7ig8tsAM8CbTPT1gDPJLWvDrPblgInw2WyjcCNZjYuTCq4EdgY1tWZ2dIwm211h32lGkNEhpg9h+soe6+GVUtnk6u6bBmVF8egZjYC+DTwH5Kavw1sMLN7gf3AytD+HHALsJdoBtyXANy9xsy+BWwL/b7p7jVh+cvAj4DhwK/Cq6sxRGSIWVeWoCA3hy+UqC5bpsWSeNy9HpjQoe0Y0Sy3jn0d+Eon+3kEeCRFezlwaYr2lGOIyNByqqGZp145wL/76FQmjCqMO5whJ+5ZbSIiGfezVw9wqqGZe1SXLRZKPCIypLg768oSLJxWxJUzVZctDko8IjKkbNt3nLcO1XHP0tmEalqSYUo8IjKklJYlGD0sjxVXTO++s6SFEo+IDBlH6s7y/M4qVi6ayfAC1WWLixKPiAwZT2ytoKnFWbVUddnipMQjIkNCc0srP966n0/Mm8gFxaPiDmdIU+IRkSHhxbeOUHXyrB72lgWUeERk0Dt6qoFH/uV9po0ZxrL5qssWt1gqF4iIpIO7s7+mnt0Ha9l1sJZdB0+yu6qWw7UNAHzj5vnk5er37bgp8YjIgNTY3MqeI3XtSWb3wVrerKqlrqEZgNwc48LiUVz7kYksmFbEpdPHsGTu+JijFlDiEZEBoO5sE28dqmPXgZNRkqmqZc/hUzS2tAIwPD+XS6aO5tYrp7NgWhELpxVx0eTRepR1llLiEZGscqTubPsZzO5wuWzfsfr29RNGFrBgWhFf+vgcFk4bw8JpRcyZMFKPNhhAlHhEJBatrU6ipj66D9N+T6aWo6ca2vvMGj+CBVOL+PxVM1g4vYgFU8cwuahQpW4GOCUeEUm7huYW9hw+dU6SebOqltONLQDk5RgXThrFJy8qZuG0IhaEV9Gw/Jgjl3RQ4hGRflV7tinpMll0qWzvkVM0tzoAIwtyuWRqEbcvmsHCaWNYMK2IeZNHUZin+zFDhRKPiPSJu3O4toHdVSfZdaC2/ab//poP7sdMHFXIwmlF3DB/UrjpP4bZ40eQo/sxQ5oSj4h0q6XV2Xfs9Ad/GxPOaI6dbmzvM2fCCC6bPoY7rp7ZPrNs0uhhMUYt2UqJR0TOcbaphXcO17XPLNt18CRvHaqjPtyPyc81Lpo8mhvmT2LhtCIWTh/D/CmjGa37MdJDSjxpsnHXIf5sww4K83MozMulMC+HgrwcCvOj5eiVS2F+DsPCe3tbXs452xV2sl3H/sNCn4LcHM36kR45Wd/ErqqT59yT2Vt9ipZwP2Z0YR6XTCviCyUz22/6z5s0moI8/fW/9J0ST5rMGDeclSUzaWhu4WxTKw3NLTQ0t0avphbqzjZztLkxam8K7aFPY3PreY/f84R17vph+Z0ntk4TYYd95OeaEl+WcXeqTp495yxm18FaDpw4095nclEhC6YW8ekFk9uTzMxxuh8j/U+JJ02iP2wb06dtW1udxpakZBQS09mmlnMSVEOKhNa+nCKhtfdvauVEfWOn2zW1+Hl99xyjZ8krRZ/8vBzyc3PIz7H25YJcIz83h7zcHPJzjYLcnPA5LId+eTlGQdv2YZuOy0PhjwxbWp33qk+xu6r2nHsyx+ubADCDuRNHcuWssaxaOrv9fszEUYUxRy5DhRJPFsrJMYbl5IZyH5m/bt7S6jR2muD6JxGeamhu7598Rtjc4u1lUNIhxwjJLEpYeTlRUooSlpGXE7UnJ7uCpMSVl5T48nNzyM8z8nM+WC4ICfCDpBlt0z5mx8950Zhtyx2TZV5O12ePZxpbePtwXfsZzO6Dtbx1qJazTdExLMjN4eIpo7lp4ZT2s5j5U4oYWaj/9SU++tcnH5KbYwwvyI3t0cDuTnOr09TSSlOz09Ta2r7c2NJKc+sHy00tUbJqamlt/xy92rYPy2GbD61v385pbv/sYbtWTje2tC83h4TccR+NLa34+Z0kdik5ebWdAbZVWK48Xk+4HUPRsDwWTCvi7iWzWTC1iIXTi/hI8SjyVY1ZsowSj2QdM2v/jZ+CuKPpmZbWDye1xubW9gT6QcKKElxjh+SXcrn5g0R4TlJsaaW5pZUWh1uvmMaCUK9sxrjhurcmA4ISj0g/yM0xctsvj4pIV3QOLiIiGaXEIyIiGaXEIyIiGRVL4jGzsWb2pJm9ZWZvmtk1ZjbezDaZ2Z7wPi70NTN70Mz2mtnrZnZV0n7WhP57zGxNUvsiM3sjbPOghTuunY0hIiKZE9cZz/8Gnnf3+cDlwJvA/cCL7j4PeDF8BrgZmBde9wEPQZREgAeAJcBi4IGkRPJQ6Nu23fLQ3tkYIiKSIRlPPGZWBFwH/BDA3Rvd/QSwAng0dHsUuDUsrwDWeqQMGGtmU4GbgE3uXuPux4FNwPKwrsjdN7u7A2s77CvVGCIikiFxnPFcAFQD/8/MXjWzfzSzkcBkd68CCO+TQv/pQEXS9pWhrav2yhTtdDHGOczsPjMrN7Py6urqvn9TERH5kDgSTx5wFfCQu18JnKbrS16p/iLO+9DeY+7+sLuXuHtJcXFxbzYVEZFuxPEHpJVApbtvCZ+fJEo8h81sqrtXhctlR5L6z0zafgZwMLRf36H9N6F9Ror+dDFGp7Zv337UzBK9+H7JJgJH+7htOmVrXJC9sSmu3lFcvTMY45rd2YqMJx53P2RmFWZ2sbu/DSwDdofXGuDb4f2ZsMmzwFfN7HGiiQQnQ+LYCPyPpAkFNwLfcPcaM6szs6XAFmA18P2kfaUao6t4+3zKY2bl7l7S1+3TJVvjguyNTXH1juLqnaEWV1wlc/4j8JiZFQDvAV8iuuy3wczuBfYDK0Pf54BbgL1AfehLSDDfAraFft9095qw/GXgR8Bw4FfhBVHCSTWGiIhkSCyJx+GTBxAAAAcySURBVN1fA1Jl0WUp+jrwlU728wjwSIr2cuDSFO3HUo0hIiKZo8oF6fVw3AF0IlvjguyNTXH1juLqnSEVl3k6HyQiIiLSgc54REQko5R4REQko5R4+pGZfTcUPn3dzJ42s7Gd9FtuZm+HIqZprxdnZivNbJeZtZpZp1MjzWxfKK76mpmVZ1FcGT1eYcweFZQ1s5ZwvF4zs2fTFEuX39/MCs3sibB+i5nNSUccfYjr98ysOun4/H6G4nrEzI6Y2c5O1ndaeDjmuK43s5NJx+svMxTXTDP7tUUFm3eZ2R+n6NO/x8zd9eqnF9HfEuWF5e8A30nRJxd4l6h0UAGwA1iQ5rguAS4m+gPbki767QMmZvB4dRtXHMcrjPtXwP1h+f5U/y3DulNpjqPb7w/8EfB/wvKdwBMZOD49iev3gL/L1L+npHGvI6qOsrOT9bcQ/YmFAUuBLVkS1/XAL2I4XlOBq8LyaOCdFP8t+/WY6YynH7n7P7l7c/hYxrkVFNosBva6+3vu3gg8TlS8NJ1xvenRH+tmlR7GlfHjFWRLQdmefP/kWJ8ElrU9CiTmuGLh7i8DNV106azwcNxxxcLdq9z9lbBcR/S0gOkduvXrMVPiSZ9/zwd/uJqss+Km2cCBfzKz7WZ2X9zBBHEdrx4VlAWGhYKyZWaWjuTUk+/f3if84nMSmJCGWHobF8Dnw6WZJ81sZor1ccjm/wevMbMdZvYrM1uY6cHDZdoriaq+JOvXYxZX5YIBy8xeAKakWPUX7v5M6PMXQDPwWKpdpGg77zntPYmrB65194NmNgnYZGZvhd/S4owrLccLuo6tF7uZFY7ZBcBLZvaGu7/bH/EFPfn+aTtGXejJmD8H1rt7g5n9IdFZ2Q1pjqsn4jhePfEKMNvdT5nZLcDPiJ4nlhFmNgp4Cviau9d2XJ1ikz4fMyWeXnL3T3W13qInoX4GWObh4mgHnRU9TWtcPdzHwfB+xMyeJrqccl6Jpx/iSsvxgq5jM7MeFZRNOmbvmdlviH5b7M/E05Pv39an0szygDGk/5JOt3F5VCmkzT8Q3ffMBmn7N3U+kn/Yu/tzZvYDM5vo7mkvHmpm+URJ5zF3/2mKLv16zHSprR+Z2XLg68Bn3b2+k27bgHlmNteiWnV3EhUvjZWZjTSz0W3LRBMlUs6+ybC4jldbQVnopKCsmY0zs8KwPBG4lqjYbX/qyfdPjvV24KVOfunJaFwd7gF8lujeQTZ4FlgdZmotJRQejjsoM5vSdm/OzBYT/Xw+1vVW/TKuET2Y8013/5tOuvXvMcv0DIrB/CIqZFoBvBZebTONpgHPJfW7hWjmyLtEl5zSHdfniH5jaQAOAxs7xkU0O2lHeO3KlrjiOF5hzAlEj0ffE97Hh/YS4B/D8seAN8IxewO4N02xfOj7A98k+gUHYBjwk/DvbytwQYaOUXdx/c/wb2kH8GtgfobiWg9UAU3h39e9wB8CfxjWG/D3Ie436GKmZ4bj+mrS8SoDPpahuD5OdNns9aSfXbek85ipZI6IiGSULrWJiEhGKfGIiEhGKfGIiEhGKfGIiEhGKfGIiEhGKfGIpImZnTrP7Z8MFRG66vMb66Kyd0/7dOhfbGbP97S/SG8p8YhkoVCnK9fd38v02O5eDVSZ2bWZHluGBiUekTQLf+39XTPbadHzju4I7TmhLMouM/uFmT1nZreHze4mqVqCmT0UipHuMrP/1sk4p8zsf5nZK2b2opkVJ61eaWZbzewdM/tE6D/HzP459H/FzD6W1P9nIQaRfqfEI5J+twFXAJcDnwK+G8rJ3AbMAS4Dfh+4Jmmba4HtSZ//wt1LgI8CnzSzj6YYZyTwirtfBfwWeCBpXZ67Lwa+ltR+BPh06H8H8GBS/3LgE73/qiLdU5FQkfT7OFGV5hbgsJn9Frg6tP/E3VuBQ2b266RtpgLVSZ+/EB5VkRfWLSAqcZKsFXgiLK8Dkos9ti1vJ0p2APnA35nZFUALcFFS/yNEpYtE+p0Sj0j6dfZQtq4e1naGqAYbZjYX+DPganc/bmY/alvXjeR6WA3hvYUP/r//E6IaeZcTXf04m9R/WIhBpN/pUptI+r0M3GFmueG+y3VExTz/hehBaTlmNpno0cdt3gQuDMtFwGngZOh3cyfj5BBVpwb4Yth/V8YAVeGM6x6ix1m3uYjsqE4ug5DOeETS72mi+zc7iM5C/ou7HzKzp4BlRD/g3yF66uPJsM0viRLRC+6+w8xeJapc/B7wr52McxpYaGbbw37u6CauHwBPmdlKourRp5PW/U6IQaTfqTq1SIzMbJRHT5ycQHQWdG1ISsOJksG14d5QT/Z1yt1H9VNcLwMr3P14f+xPJJnOeETi9QszGwsUAN9y90MA7n7GzB4geq79/kwGFC4H/o2SjqSLznhERCSjNLlAREQySolHREQySolHREQySolHREQySolHREQy6v8Djg95XES7fv0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha is  1.0\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 超参数范围\n",
    "alphas = [0.01, 0.1, 1,10,100]\n",
    "# 生成ridge实例\n",
    "ridge = RidgeCV(alphas=alphas, store_cv_values=True)\n",
    "# 模型训练\n",
    "ridge.fit(X_train, y_train)\n",
    "# 预测\n",
    "y_test_pred_ridge = ridge.predict(X_test)\n",
    "y_train_pred_ridge = ridge.predict(X_train)\n",
    "# RMSE 评估模型性能\n",
    "print('The RMSE of RidgeCV on test is %.6f'  % RMSE(y_test, y_test_pred_ridge))\n",
    "print('The RMSE of RidgeCV on train is %.6f'  % RMSE(y_train, y_train_pred_ridge))\n",
    "\n",
    "# 可视化\n",
    "# print(ridge.cv_values_)\n",
    "mse_mean = np.mean(ridge.cv_values_, axis=0)\n",
    "print(np.log10(ridge.alpha_)*np.ones(3))\n",
    "plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas), 1))\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()\n",
    "print('alpha is ', ridge.alpha_)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "岭回归模型：\n",
    "- 使用L2损失+L1正则的岭回归模型，需要选定系数$\\lambda$。系数太大，则对自变量的惩罚太强，使得模型波动太小拟合程度低；系数太小。则容易产生过拟合\n",
    "- 通过设置不同数量级的$\\lambda$，记录交叉检验的值cv_values_，选取平均值最小的作为最优$\\lambda$\n",
    "- 本例中选取的$\\lambda = 1$，测试集的RMSE略低于最小二乘线性模型拟合结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\model_selection\\_split.py:1978: FutureWarning: The default value of cv will change from 3 to 5 in version 0.22. Specify it explicitly to silence this warning.\n",
      "  warnings.warn(CV_WARNING, FutureWarning)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:471: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 458775.98326939344, tolerance: 148238.15588123395\n",
      "  tol, rng, random, positive)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:471: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 2793279.4687027633, tolerance: 148238.15588123395\n",
      "  tol, rng, random, positive)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:471: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 3122546.912490785, tolerance: 148238.15588123395\n",
      "  tol, rng, random, positive)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\model_selection\\_split.py:1978: FutureWarning: The default value of cv will change from 3 to 5 in version 0.22. Specify it explicitly to silence this warning.\n",
      "  warnings.warn(CV_WARNING, FutureWarning)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:471: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 458775.98326939344, tolerance: 148238.15588123395\n",
      "  tol, rng, random, positive)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:471: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 2793279.4687027633, tolerance: 148238.15588123395\n",
      "  tol, rng, random, positive)\n",
      "C:\\Users\\dorot\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:471: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 3122546.912490785, tolerance: 148238.15588123395\n",
      "  tol, rng, random, positive)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The RMSE of Lasso on test is 800.884918\n",
      "The RMSE of Lasso on train is 745.404118\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEGCAYAAABYV4NmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhVRdbo/+/KCAQSCAkIBBIgAUQZhIDgwCiKPWGrOLaioqit2NqTevte/XX73r6Or7MoAk5NO6Hd+toKIqA4ABIQQSWQEAIJYzCMCWRcvz9ORQ8YkgDZ2cnJ+jzPeXJ27ao6aztkpar2qS2qijHGGFPfwvwOwBhjTGiyBGOMMcYTlmCMMcZ4whKMMcYYT1iCMcYY44kIvwNoLBISEjQlJcXvMIwxpklZsWLFLlVNrO6cJRgnJSWFjIwMv8MwxpgmRUQ2He2cTZEZY4zxhCUYY4wxnrAEY4wxxhOeJhgRaSsic0QkU0TWishwERkoIktFZJWIZIjI0CPaDBGRChG5OKhskohkudekoPLBIrJGRLJF5AkREVceLyLzXf35ItLOy+s0xhjzU16PYB4H5qpqH2AAsBZ4EPirqg4E7nHHAIhIOPAAMC+oLB64FzgdGArcG5QwpgFTgDT3Gu/K7wIWqGoasMAdG2OMaUCeJRgRiQVGADMBVLVUVfcACsS6anHA1qBmU4G3gJ1BZecB81W1UFV3A/OB8SLSCYhV1SUa2LHzZeAC12YC8JJ7/1JQuTHGmAbi5W3KPYAC4AURGQCsAH4H3A7ME5GHCSS4MwBEpAvwa2AMMCSony5AXtBxvivr4t4fWQ7QUVW3AajqNhHpUF2AIjKFwAiIbt26HfeFGmOM+SkvE0wEMAiYqqrLRORxAlNVccAdqvqWiFxCYIRzDvAYcKeqVrillCrCT2kN5XWmqtOB6QDp6enH9dyC5bmFfJq1i+iIMCLDhajwMCIjwogKDyPK/Yx07yODyqJc/cjwMNf2xzqR4cIR/wyMMabJ8TLB5AP5qrrMHc8hkGDOIjCSAXgTmOHepwOvuV+sCcDPRKTc9TMqqN8k4GNXnnREedV02w4R6eRGL504fMqtXq3ctJsnFmTVe7/BSaj65BR2RDKTH5LZ0ROc/CSZ/aTPI+oF160qCw+z5GeMqZ1nCUZVt4tInoj0VtV1wFjgOwJTZyMJJIkxQJar372qrYi8CLynqv92i/x/D1rYPxe4W1ULRWS/iAwDlgFXA0+6Ou8Ck4D73c93vLrOG0f2ZMqIHpRVKGUVlZSWV1JWUUmJ+1laUUlZuVJaUUFpubpjVx5cL+hnaYUedvxjedCx67f4YNlh/f3k88srqaznZ8qFhwnxMVGkJ7cjPSWeISnt6Nsplohwu+vdGPMjr7eKmQrMFpEoIAe4lsAv+8dFJAI4hFsDORqXSO4Dlruiv6lqoXt/M/Ai0BL4wL0gkFjeEJHJwGZgYr1dUTVEJDCCiAgjJtrLTzo+FZX60wRVTSIqqzh6Ijyy/dY9h1ieW8gH32wHoFVUOIO6tSM9pR1DUuIZ2LUtMdG2E5ExzZnYI5MD0tPT1fYiO3bb9h4kI3c3y3MLWZ67m8zt+1ANjHJO7Rz7wwhncHI8iW0aYfY1xpwQEVmhqunVnrMEE2AJpn7sO1TGyk0/JpxVeXsoLa8EoEdCDOkpVdNq8aS0b2U3MxjTxFmCqQNLMN4oKa/gmy37WJ5bSEZuIRmbdrOnuAyAhNbRDEmxdRxjmrKaEoxNkhtPRUeEMzi5HYOT28HInlRWKhsKDvBlbuEPU2vVreMMTYlnYLe2tIqy/0SNaapsBOPYCMY/2/YeZHnubjJqWcdJT4knobWt4xjTmNgUWR1Ygmk89h0qY8WmHxNOdes4Q9w6TrKt4xjjK0swdWAJpvEKrOPsPWyUs/fg4es4VQnn5E5tbB3HmAZkazCmSQus48QzODm+zus4Q9y0mq3jGOMfG8E4NoJp2mpbxxmSEk96SjzpKe1sHceYemRTZHVgCSa07D1YxsrNR1/HGeKSja3jGHNiLMHUgSWY0FbTOk5iG/d9nGRbxzHmWNkajGn2qlvHyS444L4AGljHeX9NYB0nJiqcQclVCcfWcYw5XjaCcWwEY462jhMRJpzSJY6LB3XhN8OSbTrNmCA2gjGmDjrFteRXA1ryqwGdgcPXcT7L2sX/eedb8vcc5K7xfSzJGFMHlmCMOYq4lpGM7t2B0b078Idxyr3vfstzn+Rw4FA59004lTB78JoxNbIEY0wdhIUJf5twCjHRETz7yQaKSyt46OL+djOAMTWwBGNMHYkId53fhzYtInho3jqKS8t54vLTiI4I9zs0Yxol+/PLmGN0y+hU/r9f9mXetzu4/qUMikvL/Q7JmEbJEowxx+GaM7vz0MX9+Tx7F5Nmfcm+Q2V+h2RMo2MJxpjjNDG9K09dMYhVeXu44vmlFBaV+h2SMY2KJRhjTsDP+nVi+tXpZO04wKXPLWHHvkN+h2RMo2EJxpgTNLp3B166bihb9xxk4rNLyCss9jskYxoFSzDG1INhPdoz+4Zh7D1YxsRnl5C984DfIRnjO0swxtSTgV3b8vqNwyivVC59bgnfbNnrd0jG+MrTBCMibUVkjohkishaERkuIgNFZKmIrBKRDBEZ6upOEJHVQeVnBfUzSUSy3GtSUPlgEVkjItki8oS4/TtEJF5E5rv680WknZfXaUyVPifF8uZNw4mOCOPy55eyYtNuv0Myxjdej2AeB+aqah9gALAWeBD4q6oOBO5xxwALgAGu/DpgBgSSBXAvcDowFLg3KGFMA6YAae413pXfBSxQ1TTX711eXqQxwbonxPDmzWeQ0Dqaq2Yu4/PsXX6HZIwvPEswIhILjABmAqhqqaruARSIddXigK3u/AH9cWvnGFcP4DxgvqoWqupuYD4wXkQ6AbGqusS1exm4wLWZALzk3r8UVG5Mg+jStiWv3ziMbvGtuPaF5cz/boffIRnT4LwcwfQACoAXROQrEZkhIjHA7cBDIpIHPAzcXdVARH4tIpnAfwiMYgC6AHlB/ea7si7u/ZHlAB1VdRuA+9mhugBFZIqbjssoKCg4sas15ggd2rTgtSnDOLlzLDf9YwXvrNrid0jGNCgvE0wEMAiYpqqnAUUEpqpuBu5Q1a7AHbgRDoCq/stNp10A3OeKq9uyVmsorzNVna6q6aqanpiYeCxNjamTtq2imH396QxJacftr6/i1S83+x2SMQ3GywSTD+Sr6jJ3PIdAwpkEvO3K3iSwrnIYVV0M9BSRBNdP16DTSQSm1fLd+yPLAXa4KTTcz531cUHGHI/W0RG8eO1QRvVK5O631zDj0xy/QzKmQXiWYFR1O5AnIr1d0VjgOwJJYKQrGwNkAYhIatBdYIOAKOB7YB5wroi0c4v75wLz3NTXfhEZ5tpdDbzj+n2XQCLD/awqN8YXLSLDee6qdH7erxP/9Z+1PPbReuxpsibUeb1d/1RgtohEATnAtQR+2T8uIhHAIQJ3gQFcBFwtImXAQeBSt3hfKCL3Actdvb+paqF7fzPwItAS+MC9AO4H3hCRycBmYKJ3l2hM3URFhPHE5afRKiqcxz7K4sChcv7y85Pt6ZgmZIn9FRWQnp6uGRkZfodhmoHKSuVv733Hi1/kcvnQrvzXBf0It6djmiZKRFaoanp15+yBY8Y0sLAw4d5f9qV1dARPLcqmqKSCRy4ZQKQ9HdOEGEswxvhARPjjeb1p3SKC+z/IpLi0nKeuGESLSHs6pgkd9ieTMT66aWRP7rvgVD5au5PJLy2nqMSejmlChyUYY3x21bBk/vuSASzNKeSqmcvYe9CejmlCgyUYYxqBCwcl8fQVg1izZS+XT1/KrgMlfodkzAmzBGNMIzH+1JOYMWkIObsOcMlzS9i296DfIRlzQizBGNOIjOyVyCuTT6dgXwkTn13Cpu+L/A7JmONmCcaYRmZISjz/vGEYRSXlTHx2Cet37Pc7JGOOiyUYYxqhfklxvH7jcAAufW4Ja/Lt6Zim6bEEY0wj1atjG968aTitoiK44vmlLM8trL2RMY2IJRhjGrHk9jHMuXk4ibGBp2MuXm/PLTJNhyUYYxq5TnEteePG4XRPaM31L2Uw95vtfodkTJ1YgjGmCUhoHc1rNwzj1C6x3PLPlfzrq/zaGxnjM0swxjQRca0ieWXy6ZzePZ47Xv+aV5Zu8jskY2pkCcaYJiQmOoJZ1wzhnJM78H/+/Q3PfrLB75CMOSpLMMY0MS0iw5n2m8H8ckBn7v8gk4fnrbOnY5pGybbrN6YJigwP47FLBxITFc5Ti7I5UFLOPb/oS5g9uMw0IpZgjGmiwsOE/3dhP2KiI5j52UaKSsq5/6L+9nRM02hYgjGmCRMR/vfPT6Z1dASPL8iiuLSCRy8dSFSEzX4b/1mCMaaJExHuGNeL1tER/N/311JcWs603wy2p2Ma39mfOcaEiBtG9ODvv+7Hx+sLuOaFLzlgT8c0PrMEY0wIueL0bjx26UCW5+7myhnL2FNc6ndIphmzBGNMiJkwsAvP/mYwa7fu47LpSynYb0/HNP6wBGNMCBrXtyOzrhnCpu+LueS5JWzZY0/HNA3P0wQjIm1FZI6IZIrIWhEZLiIDRWSpiKwSkQwRGerqXikiq93rCxEZENTPeBFZJyLZInJXUHl3EVkmIlki8rqIRLnyaHec7c6neHmdxjRGZ6Ul8I/rh7LrQAmXPLuEjbvs6ZimYXk9gnkcmKuqfYABwFrgQeCvqjoQuMcdA2wERqpqf+A+YDqAiIQDTwPnA32By0Wkr2vzAPCoqqYBu4HJrnwysFtVU4FHXT1jmp3ByfG8esMwDpZVMPHZJWRu3+d3SKYZ8SzBiEgsMAKYCaCqpaq6B1Ag1lWLA7a681+o6m5XvhRIcu+HAtmqmqOqpcBrwAQREWAMMMfVewm4wL2f4I5x58e6+sY0O6d2ieONG4cRESZc+txSVuXt8Tsk00x4OYLpARQAL4jIVyIyQ0RigNuBh0QkD3gYuLuatpOBD9z7LkBe0Ll8V9Ye2KOq5UeUH9bGnd/r6h9GRKa4abqMggJ7kJMJXakdAk/HjG0ZwZXPL2Vpzvd+h2SaAS8TTAQwCJimqqcBRcBdwM3AHaraFbgDN8KpIiKjCSSYO6uKqulbayivqc3hBarTVTVdVdMTExNrvyJjmrCu8a1488Yz6NS2JZNmfcmidTv9DsmEOC8TTD6Qr6rL3PEcAglnEvC2K3uTwBQYACLSH5gBTFDV74P66RrUbxKBabVdQFsRiTii/LA27nwcYA80N83eSXEteOPG4aR1bM2UlzP4z+ptfodkQphnCUZVtwN5ItLbFY0FviOQBEa6sjFAFoCIdCOQeK5S1fVBXS0H0twdY1HAZcC7GtiffBFwsas3CXjHvX/XHePOL1Tbz9wYAOJjovjnDcMYkNSWqa+u5M2MvNobGXMcvN6LbCow2yWGHOBaAkngcTeyOARMcXXvIbBO8oxbjy9301flInIrMA8IB2ap6reuzZ3AayLyX8BX/DjdNhN4RUSyCYxcLvP4Oo1pUmJbRPLy5KHc+MoK/jRnNUUl5VxzZne/wzIhRuwP+4D09HTNyMjwOwxjGlRJeQVT//kVH363gz+d15tbRqf6HZJpYkRkhaqmV3fOvslvTDMWHRHOM1cO4tendeGheet4YG6mPR3T1Bvbrt+YZi4iPIxHJg6gVVQ40z7ewIFD5fz1V6fY0zHNCbMEY4whLEz4rwtOpXV0BM8tzqGotJwHL+pPRLhNcpjjZwnGGAMEHlx21/l9aB0dwSPz11NcUsHjlw8kOsIeXGaOj/15Yoz5gYgwdWwa9/yiL3O/3c49//629kbGHIUlGGPMT1x3VneuO7M7c1bmk2u7MJvjZAnGGFOtm0b2IDxMeObjbL9DMU2UJRhjTLU6xLbgiqHdeHvlFvIKi/0OxzRBlmCMMUd148gehIkw7ZMNfodimiBLMMaYo+oU15KJ6Um8mZHHVnvssjlGlmCMMTW6eVRPVOE5G8WYY2QJxhhTo6R2rbhoUBKvLs9j575DfodjmhBLMMaYWv12dE8qKpXnFuf4HYppQizBGGNqldw+hgkDOzN72SZ2HSjxOxzTRFiCMcbUyS2jUykpr+T5T20UEypUlWc/2cC2vd7cwGEJxhhTJz0TW/PL/p15ZckmCotK/Q7H1IOP1xVw/weZzP9uhyf9W4IxxtTZrWNSKS6tYNZnG/0OxZygykrlgbmZJLdvxWVDunnyGZZgjDF11qtjG37W7yRe+iKXvcVlfodjTsA7X28hc/t+/nBub6IivEkFlmCMMcfk1tFp7C8p54UvbBTTVJWUV/DIh+s5pXMsv+jXybPPsQRjjDkmfTvHMq5vR2Z9tpH9h2wU0xTNXrqZ/N0HuXN8H0+fXGoJxhhzzG4bk8a+Q+W8vGST36GYY7T/UBlPLcrmjJ7tOTstwdPPsgRjjDlm/ZLiGN07kRmf5lBUUu53OOYYPP/pRgqLSrlzfB9EvBu9wDEkGBE5S0Sude8TRaS7d2EZYxq7qWPT2F1cxj+W2iimqSjYX8KMT3P4eb9ODOja1vPPq1OCEZF7gTuBu11RJPCPOrRrKyJzRCRTRNaKyHARGSgiS0VklYhkiMhQV7ePiCwRkRIR+eMR/YwXkXUiki0idwWVdxeRZSKSJSKvi0iUK492x9nufEpdrtMYU3eDurXj7LQEnv80h4OlFX6HY+rgqYVZlJRX8odzezXI59V1BPNr4FdAEYCqbgXa1KHd48BcVe0DDADWAg8Cf1XVgcA97higELgNeDi4AxEJB54Gzgf6ApeLSF93+gHgUVVNA3YDk135ZGC3qqYCj7p6xph6NnVMGrsOlPLPLzf7HYqpxabvi5i9bDOXDulKj8TWDfKZdU0wpaqqgAKISExtDUQkFhgBzARQ1VJV3eP6iHXV4oCt7vxOVV0OHHlbylAgW1VzVLUUeA2YIIHJwzHAHFfvJeAC936CO8adHyteTzYa0wwN7R7PsB7xPPfJBg6V2SimMXvkw/VEhAu/G5vWYJ9Z1wTzhog8B7QVkRuAj4Dna2nTAygAXhCRr0RkhktMtwMPiUgegdHK3TV1AnQB8oKO811Ze2CPqpYfUX5YG3d+r6t/GBGZ4qbpMgoKCmoJwxhTndvGpLFzfwlvZOTVXtn44pste3n3661MPqs7HWNbNNjn1inBqOrDBEYCbwG9gXtU9clamkUAg4Bpqnoagem1u4CbgTtUtStwB26EU4PqRh5aQ3lNbQ4vUJ2uqumqmp6YmFhLGMaY6gzv2Z705HZM+3gDJeU2immMHpy3jratIrlxZM8G/dy6LvLHAAtV9U8ERi4tRSSylmb5QL6qLnPHcwgknEnA267sTQJTYLX10zXoOInAtNouAiOqiCPKD2vjzscRWOMxxtQzEWHq2DS27T3EWyu2+B2OOcIX2btYvL6AW0alEtuitl/b9auuU2SLgWgR6UJgeuxa4MWaGqjqdiBPRHq7orHAdwSSwEhXNgbIquWzlwNp7o6xKOAy4F23JrQIuNjVmwS8496/645x5xe6+sYYD4xIS2BA17Y883E2ZRWVfodjHNXAhpad41pw1fDkBv/8uiYYUdVi4ELgSVX9NYE7umozFZgtIquBgcDfgRuAR0Tka3c8BUBEThKRfOD3wP8WkXwRiXVrKLcC8wjchfaGqn7r+r8T+L2IZBNYY6mabpsJtHflvycwNWeM8YiIcNuYVPJ3H+RfX9koprH44JvtfJ2/l9vH9aJFZHiDf35E7VUAEBEZDlzJj7cC19pWVVcB6UcUfwYMrqbudgLTXNX18z7wfjXlOVQzxaaqh4CJtcVnjKk/Y/p04JTOsTyzKJsLT+tCRLhtFOKn8opKHp63jl4dW3PRoGp/tXqurv8F/I7AKOBtVf3WfYt/oXdhGWOaGhFh6pg0cr8v5r3V2/wOp9l7IyOfnF1F/Om8PoR7uKFlTeqaYIqBSgJfclxNYI1jtGdRGWOapHP7dqTPSW14cmEWFZW27OmXg6UVPPbRetKT23HOyR18i6OuCWY2MIvAGswvgV+4n8YY84OwMOHWMalsKCjig29sFOOXWZ9vZOf+Eu483/sNLWtS1wRToKr/o6obVXVT1cvTyIwxTdL5p3YitUNrnlyQTaWNYhrcnuJSnv1kA2P7dGBISryvsdQ1wdzrvol/uYhcWPXyNDJjTJMUHibcOjqVdTv28+F3O/wOp9l55uMNHCgp58/j+/gdSp0TzLUEbjMeT2BqrGqazBhjfuIX/TuR0r4VTy7Mwr6C1nC27jnIi1/kcuFpSfQ+qS77EXurrrcpD1DVfp5GYowJGRHhYdwyOpU/zVnNwsydjD25o98hNQuPfbQeFO4Y13AbWtakriOYpUFb5BtjTK0uOK0LSe1a8sQCG8U0hKwd+5mzIp+rhieT1K6V3+EAdU8wZwGr3EO/VovIGne7sjHGVCvSjWK+zt/L4qxdfocT8h6ct46YqAhuGZ3qdyg/qOsU2XhPozDGhKSLBiXx5IIsnliQxYi0BF9vmQ1lKzYVMv+7Hfzx3F7Ex0T5Hc4P6rpd/6bqXl4HZ4xp2qIiwrh5VE9WbNrNkg3f+x1OSFJVHvhgHYltornurO5+h3MY2yzIGOOpield6dAmmscX1LZxujkei9bt5MvcQm4bm0arqLpOSjUMSzDGGE+1iAznppE9WbaxkGU5NoqpTxWVgdFLSvtWXDaka+0NGpglGGOM5y4f2o2E1lE8uTDb71BCyr+/2sK6Hfv5w7m9iWyEu1c3voiMMSGnZVQ4U0b04LPsXazYtNvvcEJCSXkF/z1/Pf26xPHzfp38DqdalmCMMQ3iytOTadcqkicX2lpMffjH0s1s2XOQO8f3Icyn7fhrYwnGGNMgYqIjuP7sHny8roDV+Xv8DqdJ23eojKcWZnFWagJnpSX4Hc5RWYIxxjSYq4cnE9cykicW2FrMiXh+cQ67i8u4sxFsaFkTSzDGmAbTpkUk153ZnY/W7uDbrXv9DqdJ2rn/EDM+3cjP+3eiX1Kc3+HUyBKMMaZBXXNmCm2iI3jK7ig7Lk8uyKasopI/ntvb71BqZQnGGNOg4lpGcs2ZKXzwzXbWbd/vdzhNSu6uIl79cjOXDe1K94QYv8OplSUYY0yDu+7M7sREhfPUIhvFHItH5q8nMjyM28Y2ju34a2MJxhjT4NrFRHHV8BTeW72V7J0H/A6nSViTv5f/+Xork8/qToc2LfwOp04swRhjfHH92d2JjgjjGRvF1MmD8zJp1yqSKSN7+B1KnXmaYESkrYjMEZFMEVkrIsNFZKCILBWRVSKSISJDXV0RkSdEJNs9c2ZQUD+TRCTLvSYFlQ92z6bJdm3FlceLyHxXf76ItPPyOo0xxy6hdTS/OT2Zd77eSu6uIr/DadQ+z97Fp1m7uGV0KrEtIv0Op868HsE8DsxV1T7AAGAt8CDwV1UdCNzjjgHOB9LcawowDQLJArgXOB0YCtwblDCmubpV7aqeW3MXsEBV04AF7tgY08hMGdGD8DDhmY9tFHM0qsoDczPp0rYlvxmW7Hc4x8SzBCMiscAIYCaAqpaq6h5AgVhXLQ7Y6t5PAF7WgKVAWxHpBJwHzFfVQlXdDcwHxrtzsaq6RAPPY30ZuCCor5fc+5eCyo0xjUiH2BZcMbQbb6/cQl5hsd/hNErvr9nO6vy93DGuFy0iw/0O55h4OYLpARQAL4jIVyIyQ0RigNuBh0QkD3gYuNvV7wLkBbXPd2U1ledXUw7QUVW3AbifHaoLUESmuGm6jIKCguO/UmPMcbtxZA/CRHj2kw1+h9LolFVU8tC8THp3bMOvT+tSe4NGxssEEwEMAqap6mlAEYGpqpuBO1S1K3AHboQDVLdbmx5HeZ2p6nRVTVfV9MTExGNpaoypJ53iWjIxPYk3M/LZtveg3+E0Kq8vzyP3+2L+dF5vwhvphpY18TLB5AP5qrrMHc8hkHAmAW+7sjcJrKtU1Q9+Yk4SgemzmsqTqikH2OGm0HA/d9bD9RhjPHLzqJ5UqvLcJzl+h9JoFJeW8/iCLNKT2zH25GonYRo9zxKMqm4H8kSkaj+DscB3BJLASFc2Bqjau/td4Gp3N9kwYK+b3poHnCsi7dzi/rnAPHduv4gMc3ePXQ28E9RX1d1mk4LKjTGNUFK7Vlw0KIl/frmZnfsO+R1Oo/DC57kU7C/hrvP74G6QbXK8votsKjBbRFYDA4G/AzcAj4jI1+54iqv7PpADZAPPA78FUNVC4D5guXv9zZVBYLpthmuzAfjAld8PjBORLGCcOzbGNGK/Hd2Tikpl+mIbxewuKuXZjzdwzskdSU+J9zuc4xbhZeequgpIP6L4M2BwNXUVuOUo/cwCZlVTngGcWk359wRGTMaYJiK5fQwTBnbmH8s2cdOoniS0jvY7JN88vSibotJy/jy+8W9oWRP7Jr8xptG4ZXQqJeWVzPh0o9+h+GbLnoO8vGQTFw5KolfHNn6Hc0IswRhjGo2eia35Zf/OvLwkl91FpX6H44tH568HgTvG9fI7lBNmCcYY06jcOiaV4tIKZn3e/EYx63fs5+2V+UwankyXti39DueEWYIxxjQqvTq24Wf9TuLFz3PZe7DM73Aa1INz1xETFcFvR6X6HUq9sARjjGl0bh2dxv6Scl78PNfvUBrM8txCPlq7g5tG9aRdTJTf4dQLSzDGmEanb+dYzjm5IzM/y2H/odAfxagqD3yQSWKbaK49M8XvcOqNJRhjTKN029hU9h0q5+Ulm/wOxXML1u4kY9Nufjc2jVZRnn57pEFZgjHGNEr9k9oyqnciMz7Noaik3O9wPFNRqTw4L5PuCTFcOqRr7Q2aEEswxphGa+qYNHYXlzF7WeiOYv711RbW7zjAH8/tTWR4aP1KDq2rMcaElMHJ7TgrNYHpi3M4WFrhdzj17lBZBY/OX0//pDh+1u8kv8Opd5ZgjDGN2m1j09h1oJRXv9zsdyj17h9LN7Flz0HuHN90N7SsiSUYY0yjNrR7PKd3j+fZTzZwqCx0RjH7DpXx1KJszk5L4MzUBL/D8YQlGGNMo/e7sWns3F/Cmxl5tVduIrRqPDwAABFXSURBVKZ/ksOe4jLuHN/H71A8YwnGGNPoDe/ZnsHJ7Zj28QZKyyv9DueE7dx3iJmfbeSXAzpzapc4v8PxjCUYY0yjJyLcNjaNrXsP8dbKfL/DOWFPLMyirKKSP4TAhpY1sQRjjGkSRqQlMCApjqcXZVNW0XRHMRt3FfHql3lcPrQbKQkxfofjKUswxpgmoWoUk7/7IP/+aovf4Ry3hz9cR1R4GFPHhsaGljWxBGOMaTLG9OnAKZ1jeXpRNuVNcBSzJn8v/1m9jRvO7k6HNi38DsdzlmCMMU2GiDB1TBq53xfz3uptfodzzB6Ym0l8TBQ3jOjhdygNwhKMMaZJObdvR3p3bMOTC7OoqFS/w6mzT7MK+Cx7F7eMTqVNi0i/w2kQlmCMMU1KWJgwdWwqGwqK+OCbpjGKqaxUHpibSZe2LfnNsG5+h9NgLMEYY5qc80/tRM/EGJ5amE1lExjF/GfNNr7Zso/fj+tFdES43+E0GEswxpgmJzwssBaTuX0/H363w+9walRWUckjH66jz0ltuOC0Ln6H06A8TTAi0lZE5ohIpoisFZHhIvK6iKxyr1wRWeXqRonICyKyRkS+FpFRQf0MduXZIvKEuF3hRCReROaLSJb72c6Vi6uXLSKrRWSQl9dpjGl4v+jfiZT2rXhyYRaqjXcU89ryPHK/L+bP43sTHhZ6G1rWxOsRzOPAXFXtAwwA1qrqpao6UFUHAm8Bb7u6NwCoaj9gHPCIiFTFNw2YAqS513hXfhewQFXTgAXuGOD8oLpTXHtjTAiJCA/jltGpfLt1Hwszd/odTrWKSsp5/KMshqbEM7p3B7/DaXCeJRgRiQVGADMBVLVUVfcEnRfgEuBVV9SXQJJAVXcCe4B0EekExKrqEg38mfIycIFrMwF4yb1/6YjylzVgKdDW9WOMCSEXnNaFpHYteWJhdqMcxcz6bCO7DpRw5/mhuR1/bbwcwfQACoAXROQrEZkhIsH7IpwN7FDVLHf8NTBBRCJEpDswGOgKdAGCNx/Kd2UAHVV1G4D7WfUnQhcg7yhtjDEhItKNYr7O28PirF1+h3OYwqJSnlucw7i+HRmc3M7vcHzhZYKJAAYB01T1NKCIH6ewAC7nx9ELwCwCiSADeAz4AigHqkv7tf2pUqc2IjJFRDJEJKOgoKCWLo0xjdFFg5LoHNeCJxY0rrWYpxdlU1xazp/P6+13KL7xMsHkA/mquswdzyGQcBCRCOBC4PWqyqparqp3uPWZCUBbIMv1kxTUbxKw1b3fUTX15X5WTcTmExj9VNfmB6o6XVXTVTU9MTHxhC7WGOOPqIgwbh7VkxWbdrNkw/d+hwNA/u5iXlmyiYsHJ5HWsY3f4fjGswSjqtuBPBGpSt9jge/c+3OATFX9YepLRFpVTaGJyDigXFW/c1Nf+0VkmFu3uRp4xzV7F5jk3k86ovxqdzfZMGBv1VSaMSb0TEzvSoc20TyxMKv2yg3g0flZIHD7OaG9HX9tvL6LbCowW0RWAwOBv7vyyzh8egwC6ycrRWQtcCdwVdC5m4EZQDawAfjAld8PjBORLAJ3nt3vyt8Hclz954Hf1uM1GWMamRaR4dw0sidLcwr5cmOhr7Fkbt/H21/lc80ZKXRu29LXWPwmjWnO0k/p6emakZHhdxjGmON0sLSCsx9cyMmdYnll8um+xTH5xeV8mVvIp38eTdtWUb7F0VBEZIWqpld3zr7Jb4wJCS2jwpkyogefZu1i5ebdvsSwPLeQBZk7uWlkz2aRXGpjCcYYEzKuPD2Zdq0ieXJBw6/FqCr3f5BJhzbRXHdm9wb//MbIEowxJmTEREdw/dk9WLSugNX5e2pvUI8+WruTFZt2c/s5vWgZ1Xw2tKyJJRhjTEi5engycS0jeXJhdoN9ZkWl8uDcTHokxHBJelLtDZoJSzDGmJDSpkUk153Znfnf7eC7rfsa5DPfWplP1s4D/PG83kSE26/VKvZPwhgTcq45M4U20RE8tcj7tZhDZRU8Nn89A5LiOP/Ukzz/vKbEEowxJuTEtYxk0hkpvL9mO+t37Pf0s15Zsomtew812w0ta2IJxhgTkiaf1Z1WUeE85eFazN6DZTz9cTYjeiVyRs8Ezz6nqbIEY4wJSe1iorhqeDL/s3orGwoOePIZz32ygT3FZc16Q8uaWIIxxoSsG87uQXREGE8vqv9RzI59h5j1+UZ+NaAzp3aJq/f+Q4ElGGNMyEpoHc2VpyfzzqqtbPq+qF77fnxBFuUVyh/Obd4bWtbEEowxJqTdOKIH4WHCM4s21FufOQUHeH15Hlee3o3k9jG1N2imLMEYY0Jah9gWXD6kK2+tzCevsLhe+nzkw/VER4Rx65i0eukvVFmCMcaEvJtG9SRMhGc/OfFRzNd5e/jPmm1cf3YPEttE10N0ocsSjDEm5HWKa8nF6Um8mZHPtr0Hj7sfVeWBuZnEx0Rxw9m2oWVtLMEYY5qFm0f2pFKV5z7JOe4+Ps3axRcbvufW0am0aRFZj9GFJkswxphmoWt8Ky4c1IVXv9zMzn2Hjrl9ZWVg9JLUriVXDuvmQYShxxKMMabZuGV0KuWVyvTFxz6KeW/NNr7duo8/nNuL6Ajbjr8uLMEYY5qN5PYxTBjQmdnLNrPrQEmd25WWV/LwvHX0OakNEwZ08TDC0GIJxhjTrNwyJpVD5RXM+HRjndu8tnwzmwuLuXN8H8LCbEPLurIEY4xpVnomtuYX/TvzypJcdheV1lq/qKScJxZkMbR7PKN6J3ofYAixBGOMaXamjkmlqLSCWZ/XPoqZ+dlGdh0o5S7bjv+YWYIxxjQ7vTq24fxTT+LFz3PZe7DsqPW+P1DC9MU5nHdKRwZ1a9eAEYYGSzDGmGbp1jGp7C8p58XPc49a56lF2RSXlvMn247/uHiaYESkrYjMEZFMEVkrIsNF5HURWeVeuSKyytWNFJGXRGSNq3t3UD/jRWSdiGSLyF1B5d1FZJmIZLl+o1x5tDvOdudTvLxOY0zTc0rnOM45uSOzPt/I/kM/HcXkFRYze+lmJg7uSmqHNj5E2PR5PYJ5HJirqn2AAcBaVb1UVQeq6kDgLeBtV3ciEK2q/YDBwI0ikiIi4cDTwPlAX+ByEenr2jwAPKqqacBuYLIrnwzsVtVU4FFXzxhjDnPb2FT2Hizj5SWbfnLu0fnrEYHbx9mGlsfLswQjIrHACGAmgKqWquqeoPMCXAK86ooUiBGRCKAlUArsA4YC2aqao6qlwGvABNd+DDDHtX8JuMC9n+COcefHiq3OGWOO0D+pLaN6JzLzs40UlZT/UL522z7+tWoL15yZQqe4lj5G2LR5OYLpARQAL4jIVyIyQ0SCH5xwNrBDVbPc8RygCNgGbAYeVtVCoAuQF9Qu35W1B/aoavkR5QS3cef3uvqHEZEpIpIhIhkFBQUnfMHGmKZn6pg0CotKmb3sx1HMQ/PW0SY6gt+OTPUxsqbPywQTAQwCpqnqaQSSx11B5y/nx9ELBEYqFUBnoDvwBxHpAVQ38tAayqnl3I8FqtNVNV1V0xMT7f52Y5qjwcntOCs1gemLN3KwtIJlOd+zMHMnN49KJa6VbWh5IrxMMPlAvqouc8dzCCQc3DTYhcDrQfWvILBeU6aqO4HPgXTXT9egeknAVmAX0Nb1FVxOcBt3Pg4orNerM8aEjNvGprHrQAmvfrmZ++dm0jE2mmvOSPE7rCbPswSjqtuBPBGpur9vLPCde38OkKmq+UFNNgNjJCAGGAZkAsuBNHfHWBRwGfCuqiqwCLjYtZ8EvOPev+uOcecXuvrGGPMTQ7vHc3r3eB6Ym8lXm/dw+zm9aBllG1qeKK/vIpsKzBaR1cBA4O+u/DIOnx6DwJ1irYFvCCSVF1R1tVtDuRWYB6wF3lDVb12bO4Hfi0g2gTWWma58JtDelf+ew6fmjDHmJ343No2S8kp6JMYwcXCS3+GEBLE/7APS09M1IyPD7zCMMT5RVZ5YkM1Zae0ZnBzvdzhNhoisUNX06s5FVFdojDHNjYjwu3PsOy/1ybaKMcYY4wlLMMYYYzxhCcYYY4wnLMEYY4zxhCUYY4wxnrAEY4wxxhOWYIwxxnjCEowxxhhP2Df5HREpAH761KG6SSCw+WZzYtfcPNg1Nw8ncs3JqlrtdvSWYOqBiGQcbauEUGXX3DzYNTcPXl2zTZEZY4zxhCUYY4wxnrAEUz+m+x2AD+yamwe75ubBk2u2NRhjjDGesBGMMcYYT1iCMcYY4wlLMPVMRP4oIioiCX7H4jURuU9EVovIKhH5UEQ6+x2T10TkIRHJdNf9LxFp63dMXhORiSLyrYhUikjI3r4rIuNFZJ2IZItIs3jMuojMEpGdIvKNF/1bgqlHItIVGAds9juWBvKQqvZX1YHAe8A9fgfUAOYDp6pqf2A9cLfP8TSEb4ALgcV+B+IVEQkHngbOB/oCl4tIX3+jahAvAuO96twSTP16FPgz0CzunFDVfUGHMTSD61bVD1W13B0uBZL8jKchqOpaVV3ndxweGwpkq2qOqpYCrwETfI7Jc6q6GCj0qv8IrzpubkTkV8AWVf1aRPwOp8GIyP8Frgb2AqN9DqehXQe87ncQpl50AfKCjvOB032KJWRYgjkGIvIRcFI1p/4C/C/g3IaNyHs1XbOqvqOqfwH+IiJ3A7cC9zZogB6o7Zpdnb8A5cDshozNK3W55hBX3V+FIT8i95olmGOgqudUVy4i/YDuQNXoJQlYKSJDVXV7A4ZY7452zdX4J/AfQiDB1HbNIjIJ+AUwVkPki2TH8O85VOUDXYOOk4CtPsUSMizB1ANVXQN0qDoWkVwgXVVDekdWEUlT1Sx3+Csg0894GoKIjAfuBEaqarHf8Zh6sxxIE5HuwBbgMuAKf0Nq+myR35yI+0XkGxFZTWB68Hd+B9QAngLaAPPd7dnP+h2Q10Tk1yKSDwwH/iMi8/yOqb65GzduBeYBa4E3VPVbf6Pynoi8CiwBeotIvohMrtf+Q2SEb4wxppGxEYwxxhhPWIIxxhjjCUswxhhjPGEJxhhjjCcswRhjjPGEJRhjTpCIHDjB9nNEpEctdT6ubSfjutQ5on6iiMyta31jjpUlGGN8JCKnAOGqmtPQn62qBcA2ETmzoT/bNA+WYIypJxLwkPvy6RoRudSVh4nIM+6ZKu+JyPsicrFrdiXwTlAf00Qkw9X961E+54CIPCIiK0VkgYgkBp2eKCJfish6ETnb1U8RkU9d/ZUickZQ/X+7GIypd5ZgjKk/FwIDgQHAOcBDItLJlacA/YDrCXwjvsqZwIqg47+oajrQHxgpIv2r+ZwYYKWqDgI+4fD93yJUdShwe1D5TmCcq38p8ERQ/Qzg7GO/VGNqZ3uRGVN/zgJeVdUKYIeIfAIMceVvqmolsF1EFgW16QQUBB1fIiJTCPy/2YnAw69WH/E5lfz4mIB/AG8Hnat6v4JAUgOIBJ4SkYFABdArqP5OIOSfRGr8YQnGmPpztAcB1fSAoINACwC30eIfgSGqultEXqw6V4vg/Z5K3M8Kfvz/+w5gB4GRVRhwKKh+CxeDMfXOpsiMqT+LgUtFJNyti4wAvgQ+Ay5yazEdgVFBbdYCqe59LFAE7HX1zj/K54QBVWs4V7j+axIHbHMjqKuA8KBzvQg8EtmYemcjGGPqz78IrK98TWBU8WdV3S4ibwFjCfwiXw8sI/AEUAg8Q2cU8JF7GupXwLdADvD5UT6nCDhFRFa4fi6tJa5ngLdEZCKwyLWvMtrFYEy9s92UjWkAItJaVQ+ISHsCo5ozXfJpSeCX/plu7aYufR1Q1db1FNdiYIKq7q6P/owJZiMYYxrGeyLSFogC7qt60qmqHhSRewk8E35zQwbkpvH+25KL8YqNYIwxxnjCFvmNMcZ4whKMMcYYT1iCMcYY4wlLMMYYYzxhCcYYY4wn/n+uQumZcN7UvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha is: 1.0\n"
     ]
    }
   ],
   "source": [
    "alphas = [0.0001,0.001, 0.01, 0.1, 1,10]\n",
    "\n",
    "lasso = LassoCV(alphas=alphas)\n",
    "lasso.fit( X_train, y_train)\n",
    "\n",
    "lasso.fit(X_train, y_train)\n",
    "y_test_pred_lasso = lasso.predict(X_test)\n",
    "y_train_pred_lasso = lasso.predict(X_train)\n",
    "# RMSE 评估模型性能\n",
    "print('The RMSE of Lasso on test is %.6f'  % RMSE(y_test, y_test_pred_lasso))\n",
    "print('The RMSE of Lasso on train is %.6f'  % RMSE(y_train, y_train_pred_lasso))\n",
    "# # 观察预测值与真值的散点图\n",
    "# plt.close()\n",
    "# plt.figure(figsize=(4,3))\n",
    "# plt.scatter(y_train, y_train_pred_lasso)\n",
    "# plt.plot([0,1],[0,1], '--k')\n",
    "# plt.axis('tight')\n",
    "# plt.xlabel('true price')\n",
    "# plt.ylabel('predicted price')\n",
    "# plt.tight_layout()\n",
    "# plt.show()\n",
    "\n",
    "mses = np.mean(lasso.mse_path_, axis = 1)\n",
    "plt.plot(np.log10(lasso.alphas_), mses) \n",
    "#plt.plot(np.log10(lasso.alphas_)*np.ones(3), [0.3, 0.4, 1.0])\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()    \n",
    "            \n",
    "print ('alpha is:', lasso.alpha_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lasso回归模型：\n",
    "- 使用L1损失+L1正则的岭回归模型，类似岭回归模型，也需要选定系数$\\lambda$。\n",
    "- 通过设置不同数量级的$\\lambda$，记录交叉检验后的mse（mse_path_），选取平均值最小的作为最优$\\lambda$\n",
    "- 本例中选取的$\\lambda = 1$\n",
    "- 测试集的RMSE略大于最小二乘线性模型、岭回归模型拟合结果\n",
    "- 适合变量太多，需要筛选掉一些"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>columns</th>\n",
       "      <th>coef_lr</th>\n",
       "      <th>coef_ridge</th>\n",
       "      <th>coef_lasso</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>weathersit_1</td>\n",
       "      <td>1.749291e+16</td>\n",
       "      <td>807.353055</td>\n",
       "      <td>466.413537</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>weathersit_2</td>\n",
       "      <td>1.749291e+16</td>\n",
       "      <td>303.196830</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>weathersit_3</td>\n",
       "      <td>1.749291e+16</td>\n",
       "      <td>-1110.549885</td>\n",
       "      <td>-1413.106860</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>weekday_6</td>\n",
       "      <td>2.402360e+15</td>\n",
       "      <td>239.993700</td>\n",
       "      <td>18.178194</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>weekday_0</td>\n",
       "      <td>2.402360e+15</td>\n",
       "      <td>-159.664842</td>\n",
       "      <td>-373.860794</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>workingday</td>\n",
       "      <td>1.998868e+15</td>\n",
       "      <td>225.450849</td>\n",
       "      <td>2.340072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>holiday</td>\n",
       "      <td>1.998868e+15</td>\n",
       "      <td>-305.779707</td>\n",
       "      <td>-533.782955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>weekday_5</td>\n",
       "      <td>4.034916e+14</td>\n",
       "      <td>86.053817</td>\n",
       "      <td>70.818022</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>weekday_4</td>\n",
       "      <td>4.034916e+14</td>\n",
       "      <td>38.354468</td>\n",
       "      <td>26.151764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>weekday_3</td>\n",
       "      <td>4.034916e+14</td>\n",
       "      <td>2.650635</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>weekday_2</td>\n",
       "      <td>4.034916e+14</td>\n",
       "      <td>-42.135146</td>\n",
       "      <td>-36.731595</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>weekday_1</td>\n",
       "      <td>4.034916e+14</td>\n",
       "      <td>-165.252632</td>\n",
       "      <td>-157.593199</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>temp</td>\n",
       "      <td>2.785642e+03</td>\n",
       "      <td>1886.929558</td>\n",
       "      <td>2755.821947</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>yr</td>\n",
       "      <td>1.977924e+03</td>\n",
       "      <td>1973.779380</td>\n",
       "      <td>1969.577377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>atemp</td>\n",
       "      <td>1.093423e+03</td>\n",
       "      <td>1620.575825</td>\n",
       "      <td>1064.805970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>season_4</td>\n",
       "      <td>7.199798e+02</td>\n",
       "      <td>706.588321</td>\n",
       "      <td>566.770789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>season_2</td>\n",
       "      <td>6.450353e+01</td>\n",
       "      <td>91.640375</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>season_3</td>\n",
       "      <td>4.123943e+01</td>\n",
       "      <td>85.308924</td>\n",
       "      <td>-52.268361</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>season_1</td>\n",
       "      <td>-8.709105e+02</td>\n",
       "      <td>-883.537620</td>\n",
       "      <td>-1026.910307</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>hum</td>\n",
       "      <td>-1.518035e+03</td>\n",
       "      <td>-1230.338056</td>\n",
       "      <td>-1418.105814</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>windspeed</td>\n",
       "      <td>-1.534058e+03</td>\n",
       "      <td>-1378.013569</td>\n",
       "      <td>-1480.345222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>mnth_9</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>638.055264</td>\n",
       "      <td>582.797451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>mnth_5</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>299.243418</td>\n",
       "      <td>211.312209</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>mnth_3</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>198.476128</td>\n",
       "      <td>211.904773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>mnth_4</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>132.298700</td>\n",
       "      <td>87.290027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>mnth_10</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>150.600773</td>\n",
       "      <td>153.904362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>mnth_6</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>147.993020</td>\n",
       "      <td>7.303150</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>mnth_8</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>96.963291</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>mnth_2</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>-147.638207</td>\n",
       "      <td>-62.450487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>mnth_1</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>-352.051496</td>\n",
       "      <td>-247.046662</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>mnth_12</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>-430.509644</td>\n",
       "      <td>-356.734386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>mnth_11</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>-419.862647</td>\n",
       "      <td>-365.489824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>mnth_7</td>\n",
       "      <td>-1.335143e+16</td>\n",
       "      <td>-313.568599</td>\n",
       "      <td>-433.928576</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         columns       coef_lr   coef_ridge   coef_lasso\n",
       "16  weathersit_1  1.749291e+16   807.353055   466.413537\n",
       "17  weathersit_2  1.749291e+16   303.196830    -0.000000\n",
       "18  weathersit_3  1.749291e+16 -1110.549885 -1413.106860\n",
       "25     weekday_6  2.402360e+15   239.993700    18.178194\n",
       "19     weekday_0  2.402360e+15  -159.664842  -373.860794\n",
       "31    workingday  1.998868e+15   225.450849     2.340072\n",
       "30       holiday  1.998868e+15  -305.779707  -533.782955\n",
       "24     weekday_5  4.034916e+14    86.053817    70.818022\n",
       "23     weekday_4  4.034916e+14    38.354468    26.151764\n",
       "22     weekday_3  4.034916e+14     2.650635    -0.000000\n",
       "21     weekday_2  4.034916e+14   -42.135146   -36.731595\n",
       "20     weekday_1  4.034916e+14  -165.252632  -157.593199\n",
       "26          temp  2.785642e+03  1886.929558  2755.821947\n",
       "32            yr  1.977924e+03  1973.779380  1969.577377\n",
       "27         atemp  1.093423e+03  1620.575825  1064.805970\n",
       "3       season_4  7.199798e+02   706.588321   566.770789\n",
       "1       season_2  6.450353e+01    91.640375     0.000000\n",
       "2       season_3  4.123943e+01    85.308924   -52.268361\n",
       "0       season_1 -8.709105e+02  -883.537620 -1026.910307\n",
       "28           hum -1.518035e+03 -1230.338056 -1418.105814\n",
       "29     windspeed -1.534058e+03 -1378.013569 -1480.345222\n",
       "12        mnth_9 -1.335143e+16   638.055264   582.797451\n",
       "8         mnth_5 -1.335143e+16   299.243418   211.312209\n",
       "6         mnth_3 -1.335143e+16   198.476128   211.904773\n",
       "7         mnth_4 -1.335143e+16   132.298700    87.290027\n",
       "13       mnth_10 -1.335143e+16   150.600773   153.904362\n",
       "9         mnth_6 -1.335143e+16   147.993020     7.303150\n",
       "11        mnth_8 -1.335143e+16    96.963291    -0.000000\n",
       "5         mnth_2 -1.335143e+16  -147.638207   -62.450487\n",
       "4         mnth_1 -1.335143e+16  -352.051496  -247.046662\n",
       "15       mnth_12 -1.335143e+16  -430.509644  -356.734386\n",
       "14       mnth_11 -1.335143e+16  -419.862647  -365.489824\n",
       "10        mnth_7 -1.335143e+16  -313.568599  -433.928576"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性\n",
    "fs = pd.DataFrame({\"columns\":list(feat_names), \"coef_lr\":list((lr.coef_.T)), \"coef_ridge\":list((ridge.coef_.T)), \"coef_lasso\":list((lasso.coef_.T))})\n",
    "fs.sort_values(by=['coef_lr'],ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>type</th>\n",
       "      <th>RMSE_lr</th>\n",
       "      <th>RMSE_ridge</th>\n",
       "      <th>RMSE_lasso</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>train</td>\n",
       "      <td>745.108925</td>\n",
       "      <td>746.549433</td>\n",
       "      <td>745.404118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>test</td>\n",
       "      <td>798.511536</td>\n",
       "      <td>796.615615</td>\n",
       "      <td>800.884918</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    type     RMSE_lr  RMSE_ridge  RMSE_lasso\n",
       "1  train  745.108925  746.549433  745.404118\n",
       "0   test  798.511536  796.615615  800.884918"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 看看各模型的RMSE\n",
    "fs = pd.DataFrame({ \n",
    "                    \"type\":[\"test\", \"train\"]\n",
    "                   , \"RMSE_lr\": [RMSE(y_test, y_test_pred_lr), RMSE(y_train, y_train_pred_lr) ]\n",
    "                  , \"RMSE_ridge\": [RMSE(y_test, y_test_pred_ridge), RMSE(y_train, y_train_pred_ridge)]\n",
    "                  , \"RMSE_lasso\": [RMSE(y_test, y_test_pred_lasso), RMSE(y_train, y_train_pred_lasso)]\n",
    "                  })\n",
    "fs.sort_values(by=['type'],ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "各个模型的比较：\n",
    "- 最小二乘线性模型中，各个自变量系数数量级明显大于另外两种模型，有可能出现过拟合的现象。注意要在样本数量足够多的时候使用，否则容易出现过拟合\n",
    "- lasso模型中，有两个自变量系数为0，即模型会对自变量有筛选效果，剔除影响较小的自变量\n",
    "- 岭回归模型在测试数据集中的RMSE是最小的，在三个模型中有最优的表现。可以增加数据集的样本尝试进一步验证模型稳定性。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
